Last edited: March 25, 2008 

Preprint typeset using I^T^X style cmulatcapj v. 10/09/06 



THE STELLAR MASS ASSEMBLY OF GALAXIES FROM Z=0 TO Z=4. 
ANALYSIS OF A SAMPLE SELECTED IN THE REST-FRAME NEAR-INFRARED WITH SPITZER 

Pablo G. Perez-Gonzalez 1,2 , George H. Rieke 3 , Victor Villar 1 , Guillermo Barro 1 , Myra Blaylock 3 , Eiichi 
Egami 3 , Jesus Gallego 1 , Armando Gil de Paz 1 , Sergio Pascual 1 , Jaime Zamorano 1 , Jennifer L. Donley 3 

Last edited: March 25, 2008 

ABSTRACT 

Using a sample of ^28,000 sources selected at 3.6-4.5 microns with Spitzer observations of the 
HDF-N, the CDF-S, and the Lockman Hole (surveyed area: ^664 arcmin 2 ), we study the evolution 
of the stellar mass content of the Universe at 0<z<4. We calculate stellar masses and photometric 
rcdshifts, based on ^2,000 templates built with stellar population and dust emission models fitting 
the UV-to-MIR spectral energy distributions of galaxies with spectroscopic redshifts. We estimate 
stellar mass functions for different redshift intervals. We find that 50% of the local stellar mass density 
was assembled at 0<z<l (average SFR: 0.048 A^oyr^Mpc -3 ), and at least another 40% at l<z<4 
(average SFR: 0.074 «M©yr -1 Mpc~ 3 ). Our results confirm and quantify the "downsizing" scenario 
of galaxy formation. The most massive galaxies (A4>10 120 M®) assembled the bulk of their stellar 
content rapidly (in 1-2 Gyr) beyond z~3 in very intense star formation events (producing high specific 
SFRs). Galaxies with 10 n - 5 <A4<10 12 M Q assembled half of their stellar mass before z~1.5, and 
more than 90% of their mass was already in place at z^0.6. Galaxies with A4<10 115 A4q evolved 
more slowly (presenting smaller specific SFRs), assembling half of their stellar mass below z~l. About 
40% of the local stellar mass density of 10 9 '°<A^<10 11 ' M.q galaxies was assembled below z~0.4, 
most probably through accretion of small satellites producing little star formation. The cosmic stellar 
mass density at z>2.5 is dominated by optically faint (i?>25) red galaxies (Distant Red Galaxies or 
BzK sources) which account for ^30% of the global population of galaxies, but contribute at least 60% 
to the cosmic stellar mass density. Bluer galaxies (e.g., Lyman Break Galaxies) are more numerous 
but less massive, contributing less than 50% to the global stellar mass density at high redshift. 
Subject headings: galaxies: evolution — galaxies: starburst — galaxies: photometry — galaxies: 
high-redshift — infrared: galaxies 



1. INTRODUCTION 

In the last decade, our knowledge about the for- 
mation and evolution of galaxies has increased signifi- 
cantly with the advent of deep and/or wide photomet- 
ric and spectroscopic galaxy surveys carried out at dif- 
ferent wavelengths. This advance in our understand- 
ing of the evolution of the Universe is succ inctly rep- 
resented in the so-called Lilly-Madau plot (jLillv et al.1 
ll996UMadau et al.|[l996h . a dia gram showing the evolu- 
tion of the Star Formation Rate (SFR) density of the 
Universe as a function of look-back time (or redshift). 
Originally with only a few points in the diagram, it 
was clearly visible that in the last ~8 Gyr (i.e., about 
55% of its age) the Universe experienced a significant 
decrease (of about a factor of 10) in the rate at which 
new stars were created. Nowadays, there are more than 
80 da ta points in the Lilly-Madau diagram (see [Hopkins 
120041 for a nice compilation of results on this topic; 
see a lso ISchiminovich et al.l 120051 | Perez-Gonzalez et al.l 
120051 and iHopkins fc Beacoml |2006[ ) . and the picture is 
clearer at z<l, where there is just a factor of 2 scat- 
ter among the estimations coming from different surveys, 
and using different selection techniques and SFR tracers. 
At z>l, the uncertainties are larger, up to a factor of ^5, 
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but there is increasing evidence that the SFR density re- 
mained approximately constant for 4-5 Gyr (from z^l 
to z~4). 

Although the Lilly-Madau plot concentrates a large 
amount of information about the formation of struc- 
tures in the Universe, the (recent) SFR is not the best 
parameter to characterize the evolution of a galaxy, as 
it is an instantaneous parameter. Indeed, the stel- 
lar mass or the metallicity, which are closely linked to 
the star formation history, are more appropriate pa- 
rameters to follow the evolution of galaxies. Thus, 
an increasing number of studies explore the evolution 
of the cosmic comoving stellar mass density, showing 
that it has steadily increased i n the last 12 Gyr (see , 
e.g.. iBrinchmann fc Ellisl |2000L Dickinso net ali l2003bl 
Glaze brook et al.ll2004L iDrory et al.l l2005. Fon tana et al.l 
200l " see also the references given in Figure [5]) . 

Because of the increasingly large scale of cosmolog- 
ical surveys, the problem of the evolution of galax- 
ies is now being addressed by dividing the samples 
into ranges in stellar mass. In this context, the evo- 
lution of_j£alaxies seems to follow a 'downsizing' sce- 
nario (jCowie et al.lll996h , where the most massive galax- 
ies are formed first and the star formation contin- 
ues in less massive s ystems until more recent epochs 
(iHeavens et al.1 l2004t I Juneau et all 120051: iBauer et all 
20051: iPerez-Gonzalez et al.l 120051: iBundv et alj 120061 : 
Tresse et al.ll2007l) . Although the 'downsizing' picture is 
being confirmed by an increasing number of works, the 
quantification of the process is still very limited, given the 
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necessity of large samples of high redshift galaxies with 
multi-wavelength data to explore it (covering from the 
rest- frame ultraviolet to the near-infrared and beyond). 

In contrast with these observational results, classi- 
cal models of galaxy evolution assuming a Cold Dark 
Matter (CDM) Universe usually predict that the most 
massive galaxies assembled late via t he coalescence of 
small halos that form larger ones fe.g..[Kauffmann et al.l 
[1991 iBaugh et all H99l iSomerville et all 120011 ) . This 
contradicts the observational evidence of the existence 
of large galaxies at high redshifts (some of them al- 
ready harboring old stellar populations at those early 
epochs, some with significant recent star formation), 
detected by their unusually red colors (among others , 
Elston et alll988tlDev et alJll999HDickinson et al.l l200Ct 



Im et a l. 2002; iFranx et al.ll2003T) or th eir bright emission 
at sub-millimeter w avelengt hs (e.g.. ISmail et alj Il997t 



iHughes et al.l 119981 see also iBlain et al.1 120021 for a re- 
view). More recent models based on a ACDM cosmol- 
ogy succeed in predicting the early formation of mas- 
sive galaxies by introducing very large dust extinctions, 
non-standard Initial Mass Functions, and/or suppression 
of the star formation due to the quenching of cooling 
flows due to sup ernovae or Active Galactic Nuclei (e.g. , 
Cole et al.H2000b iGranato et al.ll2004 jBaugh et al.ll2005t 
Nagamine et al.1 l2005al ICroton et al.ll2006t iBower et al.l 



20MT 

In this paper, we observationally characterize the 
build-up of the stellar mass of galaxies in the last 
~12 Gyr (almost 90% of the age of the Universe) as a 
function of the stellar mass of each object. This is done 
by estimating stellar mass functions at different redshifts. 
Given that we are interested in the stellar mass assem- 
bly of galaxies, it would be convenient to analyze a sam- 
ple whose selection is based precisely on that parameter, 
the stellar mass. From studies at low and intermedi- 
ate redshift, we know that the rest-frame near-infrared 
(NIR) emission of galaxies arises mainly from relatively 
old stars that usually dominate the total stellar mass 
of galaxies, in contrast to younger stellar populations 
that may contribute little to the rest-frame NIR emission 
and stellar mass, but emit strongly at bluer wavelengths. 
Indeed, stellar mass estimations based (only) on pho- 
tometry at rest-frame wavelengths bluer than ~600 nm 
are particularly troublesome because of the ability of a 
small population of young stars to dominate the output 
of a galaxy. In the red and NIR, the light is dominated 
by similar stellar populations, but the rest-frame NIR is 
preferred for estimating stellar masses because of its rel- 
ative immunity to extinction. In addition, data at red 
wavelengths is crucial to detect galaxies that are very 
faint in the optical (too faint for optical surveys) but 
may contribute significantly or even dominate the stellar 
mass density of the Universe at high-z (e.g., Extremely 



Red Objects, EROs. lElston et al.lll988UYan et al l 



or Distant Red Galaxi es, DRGs, Franx et alj 
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Ivan Dokkum et all l2003h . These galaxies are usually 
missed by selection techniques based on rest-frame ul- 
traviolet colors (e .g., Lyman Break Galaxies, LBGs; 
iSteidel et alj 12003). Therefore, a sample selected in the 
rest-frame NIR is the most adequate to attempt a stel- 
lar mass function analysis. Still, mass-to-light ratios in 
the rest-frame NIR from galaxy to galaxy may still vary 
by factor of 6-15 (depending on the mean age of the 



stellar po pulation, the pr esence of recent burs ts, et c.; 
see, e.g.. Bell fc de Jonel 120011 IShaplev et all 120051 . or 
lLabbe et alE oOS). This means that a complete study of 
the optical-to-NIR spectral energy distribution of galax- 
ies in a galaxy-by-galaxy basis should be performed to 
obtain robust stellar mass estimates. 

This paper is based on the analysis of a sample of 
galaxies at 0<z<4 selected in 3 different fields (to min- 
imize cosmic variance problems) at 3.6-4.5 /mi with 
the Infrared Array Camera (IRAC, iFazio et al] 120 04b) 
on-bo ard of the Spitzer Space Telescope (jWerner et al.1 
2004). Even at the highest redshift in the sample, the 
sources are still selected in the rest-frame NIR (ap- 
proximately the J-band), so an IRAC selected sample 
uniquely constitutes a statistically complete sample in 
stellar mass at all redshifts up to z~4 (to a certain lower 
limit based on the flux cut of the sample). In addition, 
the estimations of the stellar masses of our galaxies al- 
ways count with a NIR band, which significantly reduces 
the uncertainties in the derived stellar masses (see, e.g., 
iFontana et al.|[2006T ) . since the relatively old stellar pop- 
ulation contributing the most to the total stellar mass 
of galaxies usually dominates the emission at NIR wave- 
lengths, and also because the NIR is relatively free of 
extinction effects and hence is better for estimating stel- 
lar masses than shorter wavelengths. Our sample selec- 
tion constitutes an extension (in area, depth, and con- 
sequently, in the number of galaxies detected) of those 
used by o ther groups based on ground-based if-ba nd 
data (e.g.- lDrorv et al.ll2004l and lFontana et al.ll2004f ). 

This paper is organized as follows: Section presents 
the dataset and samples of galaxies used in this work. 
Section [3] describes the stellar population and dust emis- 
sion models used to estimate photometric redshifts, stel- 
lar masses, and SFRs for all galaxies in our sample. Here, 
we also discuss the uncertainties in these parameters. 
Sections [4] and [5] discuss the main results about pho- 
tometric redshifts and stellar masses. More precisely, we 
present stellar mass functions and densities, discussing 
their evolution with redshift. Section [6] divides our sam- 
ple into several sub-types (such as DRGs or LBGs), and 
discusses the evolution of galaxies of different natures 
and their role on the evolution of the stellar mass den- 
sity of the Universe as a whole. Section [7] analyzes the 
SFRs of the galaxies in our sample and the evolution of 
the cosmic SFR density. Finally, Section [8] summarizes 
the conclusions of this paper. 

Throughout this paper, we use a cosmology with Ho = 
70 kms _1 Mpc- 1 , fl M = 0.3 and A = 0.7. All magni- 
tudes refer to the AB system. Th e results about stel- 
lar masses assume a ISalpete"r1 (|1955l ) universal (i.e., con- 
stant through time) Initial Mass Function (IMF) with 
0.1<.M<100 M.q and a single power-law slope in this 
range. 



2. SAMPLE SELECTION 

This paper analyzes the main properties of the galaxies 
selected by IRAC (hereafter, the IRAC selected sample), 
which should be close to a stellar mass selected sample 
up to the highest redshifts in our survey. We comple- 
mented this dataset with a sample of galaxies selected in 
a ground-based optical image (the /-band selected sam- 
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pie 4 , hereafter), in order to check the effect on our results 
of the galaxies missed by IRAC, i.e., galaxies which are 
relatively faint in the rest-frame NIR but can be detected 
in deep optical imaging. This sample of NIR-faint galax- 
ies should allow us to probe the stellar mass functions 
at small masses below the IRAC detection limits (and at 
higher masses, where the galaxies should also be detected 
by IRAC). 

The IRA C sample is drawn from th e Spitzer GTO 
see, e.g.. IPerez-Gonzalez et al.l [20051) and GOODS 
Dickinson et alj l2003ah IRACand MIPS observations 
of the Hubble Deep Field North (HDF-N) and the 
Chandra Deep Field South (CDF-S), and the Spitzer 
GTO data in the Lockman Hole Field (LHF). In each 
field, we concentrated on a relatively reduced sky area 
with the deepest coverage by Spitzer, and also observed 
by other X-ray, ultraviolet (UV), optical, near-infrared 
(NIR), and mid- infrared (MIR) surveys. In the HDF- 
N, we focused our analysis in 257 arcmin 2 centered at 
a = 12 ft 38 m 56 s , S = +62°14'06", J2000, and includ- 
ing the entire GOODS ACS footprint; in the CDF- 
S, we focused on a rectangle of 225 arcmin 2 centered 
at a = 03' i 30 m 28 s , 5 = -27°48T8", J2000, also in- 
cluding the entire GOODS ACS footprint; and in the 
LHF, we used a square area of 183 arcmin 2 centered at 
a = 10 A 52 m 47 s , 5 = +57°29'06", J2000. This adds up a 
total surveyed area of 664 arcmin 2 . 

The reduction, source extraction, and photometry of 
the IRAC and MI PS images were performed in the same 
way explained in IPerez-Gonzalez et al.l (2005) . We de- 
scribe the procedure with more details in Appendix [AJ 
The IRAC sample was built by detecting sources sepa- 
rately in the 2 bluer IRAC bands (at 3.6 /jm and 4.5 /im), 
and then merging the catalogs, and removing repeated 
sources. Aperture photometry was measured in the 4 
IRAC images (fixing the positions and forcing the detec- 
tion in all bands), obtaining the final integrated magni- 
tude after applying an aperture correction based on em- 
pirical Point Spread Functions (PSFs). All the sources in 
the IRAC sample have measured fluxes at both 3.6 /im 
and 4.5 /im. For the MIPS 24 /im images, we mea- 
sured integrated fluxes using PSF fits and aperture cor- 
rections. The /-band selected sa mple was built by detect- 
ing sources with SEXTRACTOR (jBertin fc Arnoutsl[l996h 
in the optical images. 

Our IRAC selected sample is composed of 9,074 sources 
in the HDF-N, 9,676 in the CDF-S, and 9,149 in the LHF, 
for a total of 27,899 sources (i.e., 42 sources/arcmin 2 ). 
Out of these, less than 3% (700 sources) are identified 
as stars (see the star-galaxy separation method in Sec- 
tion [A3])- Based on simulations carried out by adding 
artificial sources to the IRAC images and trying to re- 
cover their detection and input flux, we estimate that 
our IRAC catalogs in the HDF-N and the CDF-S are 75% 
(90%) complete down to 1.6 /xJy (5.0 fj,Jy) at 3.6 /im, and 
1.4 juJy (4.0 /iJy) at 4.5 /im. For the LHF, where deep 
GOODS IRAC data are not available, the 75% (90%) 
completeness levels are 2.2 fiJy (5.8 /iJy) at 3.6 /im, 
and 2.0 /zJy (4.8 /iJy) at 4.5 /im. Above the 75% com- 

4 For this selection, we chose the deepest ground-based images 
in a band common (or similar) to the 3 fields, namely, the Subaru 
/-band images in the LHF and the HDF-N, and the Subaru 7VB816 
image (close to an /-band image, and also very deep) in the CDF-S. 



pleteness flux limits, our sample has 7,512 galaxies (after 
removal of stars) in the HDF-N, 6,546 galaxies in the 
CDF-S, and 5,341 galaxies in the LHF, adding a total 
of 19,399 galaxies (29.2 sources/arcmin 2 ). Out of these, 
6,686 (35%) galaxies are detected by MIPS at 24 /im, 
3,483 (18%) above our 75% 24 /im completeness level 
[F(24)=80 /iJy]. 

We concentrated our analysis of the /-band selected 
sample on the region covered by the other UV-to-MIR 
surveys, and enclosing a similar number of sources as 
those detected with the IRAC selection (therefore, we 
considered smaller regions in each field). We used an 
area of 101 arcmin 2 centered at a = 12 ft 37 m 00 s , S = 
+62°13'30" (J2000) in the HDF-N, 103 arcmin 2 at a = 
03' l 32 m 28 s S = -27°48'54" (J2000) in the CDF-S, and 
70 arcmin 2 at a = 10 ft 52 ro 48 s 5 = +57°29'24" in the 
LHF. The samples arc formed by 7,326 sources (112 of 
them identified as stars) in the HDF-N, 6,680 (87 stars) 
in the CDF-S, and 6,797 (99 stars) in the LHF, for a total 
of 20,505 galaxies with /<25.5 (75 sources/arcmin 2 ). 

The Spitzer data were complemented with other pub- 
licly available and proprietary photometric and spectro- 
scopic data in the 3 fields. These data cover the elec- 
tromagnetic spectrum from UV to MIR wavelengths. 
The description of the different datasets and the proce- 
dure used to get merged UV-to-MIR photometry for each 
source is described in detail in Appendix lAl The spectral 
energy distributions (SEDs) of each source were used to 
remove stars from our sample, detect candidates to har- 
bor an Active Galactic Nuclei (AGN), and to estimate 
photometric redshifts, stellar masses, and SFRs for the 
entire sample in a galaxy- by-galaxy basis, as explained 
in Section [3] and Appendix [Bl 

3. ESTIMATION OF PHOTOMETRIC REDSHIFTS, 
STELLAR MASSES, AND STAR FORMATION RATES 

The estimation of the photometric redshift, stellar 
mass, and SFR of each galaxy in our IRAC and /-band 
selected samples was carried out in a two step process. 
Given the significant degeneracies inherent to any stel- 
lar population modelling, and in order to get the best 
estimations of the interesting parameters, we decided to 
first build a reference set of stellar population and dust 
emission templates. This trained template set was used 
in the second step to obtain photometric redshifts, stellar 
masses, and SFRs for the entire sample. The reference 
template set was built with the ^2,000 galaxies in our 
spectroscopic sample with highly reliable redshifts and 
well-covered SEDs (with enough data points from the 
rest-frame UV to NIR/MI R wavelengths). This is th e 
same approach we chose in iPerez- Gonzalez et"aU (|2005ft . 
As a major improvement of our photometric redshift 
technique described in that paper, we (significantly) in- 
creased the spectral resolution of the templates by fitting 
the SEDs of the galaxies in the reference spectroscopic 
sample with models of the stellar population and dust 
emission (probing more than 10 11 different models). 

In Appendix [HI we describe the stellar and dust emis- 
sion modelling procedure, the building of the reference 
template set, and the procedure to get photometric red- 
shift, stellar mass, and SFR estimates for each galaxy in 
our entire sample. In this Appendix, we also evaluate 
the goodness of our photometric redshift, stellar mass, 
and SFR estimates. We show that our photometric red- 
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shifts for galaxies at z<1.5 are better than a z / (l+z)<0.1 
(where rj z is the absolute value of <5z=z S pcc-z p hoto) for 
approximately 87% of the galaxies in our complete sam- 
ple, and better than rj z /(l+z)<0.2 for 95%. At z>1.5, 
we test our photometric redshifts distributions for dif- 
ferent samples of high redshift galaxies (LBGs, DRGs, 
and BzK sources, see Section [5] for more details), ob- 
taining acceptable results, in good agreement with other 
spectroscopic and photometric redshift analysis. 

The distributions of photometric redshift uncertainties 
(as derived from the comparison with spectroscopic red- 
shifts in Appendix [B]) for different magnitude and red- 
shift intervals are used in Section [5] to estimate the un- 
certainties in the stellar mass functions. In addition, the 
redshift intervals in that Section and the following are 
constructed assuming that the typical photometric red- 
shift error is <r z /(l+z)~0.1 (equation valid for more than 
85% of our sample). We would like to stress that the 
results described in the following Sections are more ro- 
bust at z<1.5, where the photometric redshifts are well 
tested, and the photometry is very accurate, than at 
z>1.5, where the unavailability of spectroscopic redshifts 
does not allow a characterization of the photometric red- 
shifts as thorough as at low z, and photometric errors are 
generally larger. 

In Appendix [B] we also discuss the goodness of our 
stellar mass estimates. We conclude that the choices of a 
single population or a two component population in the 
stellar emission models, the use of distinct stellar pop- 
ulation libraries, different IMFs, or different extinction 
recipes produce changes in the derived stellar masses of 
a factor of 2-3, which is also the typical error in any 
stellar population synthesis analysis linked to the degen- 
eracies of the solutions to the problem. Therefore, our 
stellar mass estimates are good within a factor of 2-3. 

The estimations of the SFRs for each galaxy are also 
proved to be good within a factor of 2 in Appendix [Bj 
an uncertainty which is consistent with other evaluations 
of UV- and IR-based SFRs Ce.g.. iPapovich fc Belll[200l 
ILe Floc'h etlH[2005l : ICatmti et al.H200l ~ 

Finally, Appendix [Bj also discusses the validity of our 
estimated parameters for galaxies harboring an AGN. We 
conclude that photometric redshifts and stellar masses 
should not be affected dramatically for most AGNs (ex- 
cept for very bright Type 1 AGNs), but IR-based SFRs 
can be overestimated. For this reason, we exclude AGNs 
from the analysis of SFRs performed in Section 17.21 but 
we keep most of them (only excluding very bright X-ray 
sources) in our calculations of the stellar mass functions 
and densities. 

4. REDSHIFT DISTRIBUTION OF OUR SAMPLE 

Figure [T] shows the photometric redshift distribution of 
our IRAC selected sample (average number density and 
number densities in each field), the subsample detected 
also by MIPS at 24 /itm, and the /-band selected sam- 
ple. Only sources with fluxes above our 75% complete- 
ness levels are included. The distributions have been 
constructed taking into account the typical photometric 
redshift error (<7 z /(l+z)~0.1), i.e, Figure[T]represents the 
real redshift distribution convolved with the photometric 
redshift probability distribution. 

Figure [T] demonstrates the importance of cosmic vari- 
ance effects on deep photometric surveys. Indeed, large 
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Fig. 1.— Redshift distribution of our IRAC, MIPS, and /-band 
selected samples (including all galaxies above the completeness 
level). For the iRAC and MIPS samples, the three fields used 
in this paper are plotted with different colors, and the average 
number densities are plotted in black. Continuous lines refer to 
the number densities for the entire IRAC selected sample (scale 
on the left vertical axis). The dashed lines refer to the subsample 
(within the IRAC sample) detected also at 24 /im (scale on the 
right vertical axis). The dash-dotted line show the redshift distri- 
bution of the /-band selected sample, with the same scale as the 
IRAC distribution. 



scale structures are clearly visible and located at different 
redshifts for our 3 fields, especially at z<l. Number den- 
sity variations of up to a factor of 2 can exist at a given 
redshift from one field to another. The HDF-N shows two 
very prominent density peaks at z^0.5 and z~0.9, con- 
sistent with t he spectroscopic re dshift histogram found 
in figure 16 of IWirth et al.l (|2004f ). There are also minor 
prominences at z=l. 5-2.0 and z=2.0-2.5, which are also 
seen i n the spectros c opic fo llow-up of UV-selected galax- 
ies in iReddv et~aT1 (|2006af) . The CDF-S presents very 
prominent density peaks at z~ 0.3, z~0.7, and z~l.l 
(the latter broadens up to z~1.4), which coincides (after 
convolution with the typical photometric redshift uncer- 
tainty) with the most prominent spectroscopica l ly con - 
firmed peaks found in figure 7 of lVanzella et all (|2006l) . 
The LHF shows an enhanced density at z<0.3, z~0. 7-1.0 
and a very prominent peak at z=1.5— 1.8. These peaks are 
co nsistent with the high d ensit y of X-ray sources found 
by iMainieri et~aTl (|2002b| ) and IZappacosta et all (|2005l ) 
at z~0.8 and z~1.6-1. 8, and the analysis of the sha llower 
IRAC SWIRE data in iRowan-Robinson et ail (|2005f ). 

Only by combining data for several fields are we able 
to smooth out cosmic variance effects. Indeed, the aver- 
age redshift distributions for IRAC sources (black con- 
tinuous line) and MIPS sources (black dashed line) are 
much smoother than the analogous curves for the in- 
dividual fields. The shape of the redshift distribution 
for the IRAC sample is typical of a flux limited sample 
with a roughly homogeneous detection probability, i.e, 
the detect ion of a source depends only on its magnitude 
(see, e.g., iBemtezI |2000|) . The detection probability of 
our IRAC survey peaks at around z=0. 8-1.0. For z<0.6, 
the detection of sources is dominated by the surveyed vol- 
ume, and after z~1.0, the detection probability decreases 
exponentially up to z^4. About half of our sample lies 
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at z>0.9, ~40% at z>l, and ~20% at z>1.5. The bulk 
of the galaxies in this study (~90%) lie at z<2. This 
implies that our results about stellar mass functions and 
densities are very robust up to z~2, just where our photo- 
metric redshifts are empirically well tested. Beyond that 
point, we still include ~3000 galaxies, enough to still ob- 
tain statistically meaningful results (although systematic 
errors such as redshift outliers will also contribute more 
to the errors above z=2). 

The statistics for the /-band selected sample are very 
similar to those for the IRAC sample: the average distri- 
bution peaks at around z=0.7, and then decays exponen- 
tially, enclosing about 50% of the sources below z=0.9, 
80% at z<1.5, and 10% at z>2.0. Figure [1] shows that 
most of the galaxies included in the /-band selected sam- 
ple and missed by IRAC lie at z<1.5. At higher redshifts, 
the number densities of the /-band and IRAC samples 
are almost identical, which is consistent with the fact 
that more than 90% of the IRAC sources were detected 
in our deep Subaru /-band images (see Appendix [Aj . 
This means that the /-band mass completeness level is 
very similar to the IRAC level, except for z<1.5, where 
the /-band should help to probe (slightly) lower masses 
than the IRAC selection (see Figure H]). We would need 
optical images deeper than /~25.5 to detect less massive 
systems at high redshift. 

Figure [1] also shows the redshift distribution of the 
IRAC sources detected by MIPS at 24 ^m and hav- 
ing F(24)=80 fiJy (dashed lines). The redshift distribu- 
tion is similar to that presented in lPerez~G onzalc z et al.l 
( 2005) , but the improvement in the photometric redshift 
estimations reveals a more pronounced density bump at 
z^l.7 and a weak bump at z~2.6. The origin of these 
bumps can be found in the increase of the detection prob- 
ability induced by prominent PAH features entering the 
MIPS 24 ixm filter as w e move to higher redshifts (see 
also iCaputi et alj 12006). Indeed, a typical PAH spec- 
trum shows an absence of features around A=10 /Ltm, 
which produces the detection local minimum at z^l.3 
observed in Figure [TJ At 6<A<10 /im, there are sev- 
eral PAH features (the most prominent at 5.5 /zm and 
7.7 /itm) that are responsible for the bumps in the red- 
shift distribution. Note that the final detected density 
for MIPS sources is a convolution of the real redshift dis- 
tribution of galaxies (affected by large scale structure), 
the detection probability (dependent on the limiting flux 
of the survey and the spectra of the galaxies), and the 
photometric redshift uncertainty distribution. These two 
effects (detection probability and redshift uncertainties) 
result in blurring out redshift-dependent features so they 
are at lower contrast to the overall real distribution. 

5. STELLAR MASS FUNCTIONS AND DENSITIES 

5.1. Completeness of the sample 

Figure [5] shows the distribution of stellar masses of in- 
dividual galaxies in our IRAC survey as a function of 
redshift. The blue continuous line shows the stellar mass 
corresponding to a passively evolving galaxy formed in a 
single instantaneous burst of star formation occurred at 
z=oo and having a 3.6 /im flux equal to the 75% com- 
pleteness level of our IRAC sample. The stellar mass cal- 
culated in this way assumes the maximum mass-to-light 
ratio given by the oldest instantaneously formed stellar 
population possible at each redshift. Any burst occur- 
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Fig. 2. — Distribution of the stellar masses of all individual 
galaxies in the IRAC (all symbols) and MIPS (red symbols) se- 
lected samples as a function of redshift (shown with a logarithmic 
scale in the quantity 1+z in the bottom horizontal axis and the 
corresponding look-back times in the top axis). The continuous 
blue line shows the stellar mass value at each redshift above which 
our IRAC survey is 75% complete for passively evolving galaxies. 
The dashed blue line shows the completeness for highly extincted 
[E(B — V)=l.l] bursts. Sources whose stellar mass is beyond the 
vertical axis scale are plotted with arrows at the source redshift. 

ring after the primary placed at z—oo should decrease the 
observed mass-to-light ratio (unless it presents a high at- 
tenuation, see below), thus giving a smaller stellar mass. 
Therefore, the values given by the blue continuous curve 
in Figure [2] are the minimum stellar masses that a max- 
imally old galaxy with a flux equal to the 3.6 ^m 75% 
completeness level should present, and our survey must 
be complete (actually, at least 75% complete) against 
passively evolving galaxies with masses above th e con- 
tinuous curve. As noted by I Font ana et al.l (|2006r i. high 
mass-to-light ratios can also be found in galaxies with 
very extincted bursts. The dashed blue line in Figure [2] 
shows the completeness level of our survey for instan- 
taneous star-forming b ursts extincted by E(B — V)=l.\ 
magnitudes (as used bv lFontana et al.l2006l based on the 
typical extinction o f highly obscur e d high redshift galax- 
ies) and following a lCalzetti et alJ (|2000l ) extinction law. 

Note that if the density of galaxies of a given stellar 
mass at a certain redshift is very small, our surveyed 
volume may not be large enough to enclose any galaxy 
of that mass (we would not detect any, although there 
might exist galaxies of that mass in the Universe at that 
redshift). This is the effect seen in Figure [2] at high stellar 
masses: at z<0.2, our surveyed volume is not enough to 
detect galaxies with AI^IO 110 A4q, and we can only 
detect galaxies with M>10 12 M® at z>0.6. It is also 
interesting to notice that the most massive galaxies with 
_A/f>10 12 M.q are only found in the regions presenting 
the highest densities, just where the redshift distribution 
for individual fields peak (see Figure [IJ. 

The estimations of the stellar mass functions in the 
following sections will be carried out for stellar masses 
above the completeness level (against passively evolving 
galaxies) shown in Figure [2] (with a continuous line), i.e., 
no completeness correction will be carried out to try to 
recover the stellar mass function at smaller masses (below 
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the blue continuous curve in Figure [2]). 

5.2. Stellar mass function estimation procedure 

The entire redshift range 0<z<4 was divided into 12 
intervals, the size of each bin chosen to have a statisti- 
cally representative number of galaxies and taking into 
account the typical photometric redshift errors. Our goal 
was to estimate stellar mass functions at each redshift 
bin. Classical met hods to obtain luminosity functions or 
mass function (see iWillmerl 1997. for a discussion about 
them) rely heavily on the use of a flux band on which 
the selection of the studied sample is based. If the band 
where the selection is based is far from the band where 
we want to estimate the luminosity function (or, in the 
case of estimating stellar mass functions, the magnitude 
is not directly and easily linked to the stellar mass of each 
gala xy), significant systematic errors a re introduced (see, 
e.g.. lLovedavll2000t lllbert et al.ll200l . In our case, our 
selection is carried out in luminosity at 3.6-4.5 //m, but 
we want to obtain a stellar mass function, which is linked 
to that luminosity (but not directly proportional). To 
solve this problem, we estimated a bivariate luminosity- 
stellar mass function for e ach redshift bin. The p r ocedur e 
is identical to that used in lPerez- Gonzale z et al.l ((2003a), 
and accounts for the fact that the selection of the sam- 
ple is carried out in a certain photometric band, while 
we eventually want the number density function relative 
to a different parameter (in our case, the stellar mass). 
The bivariate luminosity-stellar mass function, BLMF or 
&(L, Ai), is defined as the number density of galaxies (in 
a limited co-moving volume given by our surveyed area 
and the redshift interval considered) with a given lumi- 
nosity in a certain band and a given stellar mass. This 
defini tion is an exten sion of the bivariate luminosity func- 
tion (|Lovedavi [20001. The estimation of the BLMF was 
performed with a stepwise maxim um likel ihood (SWML ) 
technique (jEfstathiou et al.lll988l see also lWillmei1ll997f ). 
extended to consider two independent variables. 

To estimate stellar mass functions, we used the IRAC 
3.6 /im band as the luminosity variable in the BLMF, 
given that this is the filter where the selection of the sam- 
ple was carried out. For the /-band selected sample, we 
used the / filter as the selection band. Once the BLMF 
is estimated, if we integrate it through all luminosities, 
we can estimate the number density of galaxies with a 
given stellar mass, i.e., the stellar mass function, SMF or 
0Sm(A4). We only estimated the stellar mass function 
down to the completeness threshold of the stellar mass 
discussed in Section |5~T1 

In the classical SWML method, the errors in the BLMF 
are estimated from the covariance matrix. In our case, 
the estimation of the BLMF uncertainties was carried out 
by combining the SWML technique with a Montecarlo 
method to take into account the photometric redshift 
errors and outliers, as we did in iPerez-Gonzalez et al. 
(2005) following the procedure described in IChen et al. 
(2003). We considered the photometric redshift as a sta- 
tistical variable whose error comes from the comparison 
with spectroscopic redshifts. These errors depend on the 
actual redshift of the galaxy, so we considered different 
photometric redshift uncertainty distributions for differ- 
ent redshift intervals. We also considered the dependence 
of the redshift uncertainties on the apparent brightness of 
the source (more accurate photometry allows better es- 



timations of the photometric redshift), dividing the red- 
shift dependent distribution of redshift uncertainties into 
magnitude bins. For z>1.5, where very few spectroscopic 
redshifts are available to test our photometric redshifts, 
we only considered one single redshift and magnitude 
interval. The Montecarlo method uses the redshift un- 
certainties based on the comparison with spectroscopy, 
given that they are more reliable (they directly test the 
goodness of the photometric redshifts) than the errors de- 
rived from the probability distribution based on the x 2 
minimization, and they include the effect of outliers. The 
Montecarlo extension to the SWML method consists in 
calculating the stellar mass function by randomly varying 
the redshifts of the whole sample according to the distri- 
bution of uncertainties (which are usually non-Gaussian), 
and calculate the stellar mass function again. After 100 
iterations, the average and standard deviations of each 
point in the stellar mass function are taken as the final 
results. 

The results for the SMF (data points and uncer- 
tainties 5 ) were fitted with a smooth function using a 
iSchechterl (|1976f ) parametrization, to facilitate compar- 
ison with similar fits in the literature. Both the IRAC 
and /-band SMF estimations were used in the fits, down 
to the IRAC completeness level. For the five bins at 
highest redshifts, the faint-end slope of the SMF was 
poorly constrained by our data, so we combined our re- 
sults with other estimations of the stellar mass functions 
found in the literature. These estimations are based on 
the analysis of galaxy samples typically selected at op- 
tical wavelengths, which is more effective at probing the 
low mass regime of the stellar mass function. We only 
used literature data points at masses M>10 8 Mq at 
z<1.3 and A4>10 9 Mq at z>1.3, where the complete- 
ness levels of these optically based samples are supposed 
to be high (based on the sharp turnovers of these SMFs). 
Note that the completeness of the optically selected sam- 
ples at masses around A4~10 8-9 M.q is difficult to es- 
timate (due, for example, to the significant effect of the 
extinction, and the need of extremely deep data to probe 
this mass regime) and is not well understood (it is not 
discussed in the reference papers) , so the low mass slopes 
at z>2 should be taken with caution. The SMF points, 
errors, and Schechter fits for each redshift bin are shown 
in Figures and 2] (black filled and open stars for the 
IRAC and /-band selected samples, respectively). These 
figures also show other SMF estimations found in the lit- 
erature (color points, see captions for references). The 
plots also depict the SMFs and fits for the subsample of 
galaxies detected simultaneously by IRAC and MIPS at 
24 /im (filled circles). The data points and Schechter fit 
parameters are given in the electronic Tables Q] and 
respectively. 

5.3. The local stellar mass function and density 

Figure [3] shows our estimations of the local stellar mass 
function (including sources at 0.0<z<0.2) based on both 
the IRAC (filled stars) and /-band (open stars) selected 
samples. Given that we are surveying a very limited 

5 Note that the data point at X=10 12 Mq in each SMF, 
which accou nts for the very few high-mass sources discussed in 
Section 15. II presents a very large uncertainty (as it includes very 
few sources) and have a negligible effect on the Schechter fits. 
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Fig. 3. — Local stellar mass function estimated with the IRAC 
selected (filled stars), /-band selected (open stars), and MIPS se- 
lected (filled circles) samples at z<0.2. For clarity, the /-band data 
points have been artificially drifted from the original x-position 
(the same as the ones for the IRAC selected sample) and we do 
not show the uncertainties for the MIPS data points. The vertical 
gray dashed line shows the completeness level of our IRAC survey 
in the local Universe. The Schechter fit to the IRAC and /-band 
data (for masses .M>10 9 Mq) is shown with a black continuous 
line. Our estimation of the local stellar mass function is compared 
with th e one estimated b y lCole et all l|2001l . red crosses and line), 
and by IBell et all (120031 . blue line). The best Schechter fit to the 
data for the MIPS sample (i.e., for local star-forming galaxies) is 
plotted with a dashed line. This SMF is compared with the one 
published by Perez-Gonzalez et al. (2003a, green asterisks and line) 
for Ha-selected local star-forming galaxies. 



volume in the local Universe, we do not detect many 
sources with Ai>10 110 Mq (this explains the large er- 
rors in this mass regime), but our statistics are much 
better at low masses. Our Schechter fit only refers to the 
data points above A4>10 9 Mq, where our data is di- 
rectly comparable with previously published mass func- 
tions. Our results are very similar to tho se published by 
iCole et al.1 (f200l and IBell et all ([2001 based on NIR 
2MASS data down to the completeness limit of their sur- 
veys (A^^IO 9 5 Mq). Our deeper data co nfirm the faint- 
end slope estimated bv ICole et al.l (j2001) down to even 
smaller masses, A4~10 90 Mq. We also find a steepening 
of the stellar mass function at A4<10 9 Mq (at least for 
_A/(>10 7 - 9 Mq, our completeness level at z~0). 

By integrating our local SMF, we obtain a value for the 
local stellar mass density of p» = 10 8 - 75±0 12 Mq Mpc~ 3 , 

in excellent ag re ement w ith the v alues found in 

Salucci fc PersicT (Il999l). ICole et all <f200lh . and 
Bell et al.1 (120031) (^=10 8 - 75 ' 8 - 76 - 8 - 74 A4 Mpc- 3 , 
respectively). The good agreement of our local stellar 
mass function and density with previous estimations 
found in the literature demonstrates that aperture 
effects in our photometric catalogs are not critical, i.e., 
they do not affect our results even at low redshifts where 
the galaxies present relatively large angular sizes. The 
steepening of the stellar mass function at .M<10 90 Mq 
has no significant effect on the integrated stellar mass 



density (justifying the exclusion of these points in the 
Schechter fit): the galaxies with 10 78 <X<10 9 Mq 
contribute less than 2% to the total stellar mass density. 

Figure [3] also shows the SMF of the sources detected 
by MIPS at 24 /mi, i.e., the galaxies with active star 
formation (filled circles and dashed line fit). The MIPS 
results (data points an d fit) are in excellent agreem ent 
with those published bv lPerez-Gonzalez et alj (|2003al) for 
a Ha-selected sample of star-forming galaxies in the lo- 
cal Universe. The local stellar mass density locked in 
star-forming galaxies is p^ F = io 7 - 85±0 - 07 Mq Mpc -3 , i.e., 
13±4% of the global stellar mass density in the local 
Universe is found in active star- forming galaxies. Fig- 
ure [3] also shows that approximately 1 of every 4 galax- 
ies in the local Universe with M<10 105 Mq is forming 
stars currently and would be detected in the IR or with 
a SFR tracer such as the Ha emission. At higher stellar 
masses, the fraction of star-forming galaxies decreases by 
more than a factor of 2 (e.g., ~10% of all galaxies with 
are forming stars actively). 

5.4. The evolution of the stellar mass function 

Figure [4] presents the global stellar mass functions es- 
timated in the 12 redshift intervals up to z=4. We show 
the results obtained with the IRAC selected (filled black 
stars), the /-band selected (open black stars), and the 
MIPS selected (filled circles) samples. Other estima- 
tions found in the literature are also plotted at each red- 
shift interval (normalized to the lSalpeteilll955l IMF). Wc 
have fitted our SMF data points, and for z>1.5, also the 
data points from other works below our completeness 
level to better constrain the slope at lighter masses, to 
a Schechter function. In the case of the 24 /im galaxies, 
we assumed the same faint-end slope estimated for the 
global SMF (based on the IRAC and /-band samples). 

Our estimation of the stellar mass function is con- 
sistent (within errors and for the same mass ranges) 
with previous estimations found in the literature. It 
is i nteresting to notice that the faint-end slope derived 
by iDrorv et all {2005) at z>1.5 i s significantly larger 
than that found by other authors (jFontana et al.l [20041 : 
IConselice et al.ll2~005D . Our results are in good agreement 
with Drory's in the stellar mass range probed by both 
surveys. The possible explanations for the discrepancy 
at low masses are field-to-fiel d variations, an over estima- 
tion of the faint-end slope in IDrorv et al.l (120051) due to 
the use of the V m& ^ method (see lllbert et al.ll2004l ). or sys- 
tematic errors in the determination of the stellar masses. 
Due to this discrepancy, we did not use these points in 
our SMF fits. At z>2, we estimate number densities 
of massive galaxies whi ch can be up t o 0.8d ex higher 
than those estimated by iFontana et al.l ([2006) . Some of 
this discrepancy (up to ~20%) may be due to the under- 
density observed in CDF-S (see Section 153)) . 

Our results show that the local density of galaxies 
(shown with a gray line in all panels) with masses 
7W>10 12 M® was already reached by the SMF at z=2.5- 
3.0, i.e., the most massive galaxies were already in place 
at that redshift (approximately 11 Gyr ago). The mass 
assembly of galaxies shifts to smaller masses as we move 
to lower redshifts. By z~l, the SMF has reached nearly 
the local density for galaxies with A^IO 118 Mq. At 
z<l, the star formation in the Universe occurs mainly in 
galaxies with _M<10 n ' 5 Mq. It is also interesting to no- 
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Fig. 4. — Stellar mass functions for 12 redshift intervals from z=0 to z=4. Our estimations at each redshift interval are plotted with 
black filled stars and errors for the IRAC selected sample, and with open black stars for the /-band selected sample (errors for this sample 
are not plotted for clarit y). Filled c ircles show the SMF for galaxies detected by MIPS at 24 fim. The SMF data (our estimations and 
others) are fitted with a ISchechterl (1976) func tion (black cont inuous line for the global SMF, and dashed line for the SMF for 24 ^m 
sources). All panels show the local SMF from [Cole et al.l ( 2001) with a gray curve. The vertical dotted line shows our 75% completeness 
limit for the IRA C selected sample (continuous cur ve in Figure^. Co lor points show estimations f rom other papers: red crosses come from 
I Cole et al.l (I200H C 01); oran ge open triangles fro m lBorch et alJl2006l . B06); mag enta crosses fromlPannella et alj J2006I): red s quare s from 
iDrory et al.l H2004I D04) and lDrorv et al.l ||2005| . D05); g reen crossed circ l es fro m lFontana et al.l 12003 . FC03). IFontana et al.l | |2003. F04), 
and lFontana et al.l (120061 . F06); and p urple squares fromlCo nselice et al.l l|2005l . C05). Green asterisks at 0.0<z<0.2 show the stellar mass 
function of local star-forming galaxies ( Perez- Gonzalez et al. 2003a, P03). 



tice the significant evolution (approximately a factor of 
0.2dex or 60%, as shown by our data, a nd al so confirmed 
by th e results of iPannella et al.l 120061 and iBorch et al.l 
I2006D of the SMF between z~0.4 and z=0 (i.e., a period 
of 4 Gyr) for stellar masses 10 9 < M<10 n Mq. We w ill 
comment more on this recent evolution in Section [5.51 

Figure S] also shows that the slope of the SMF at low 
masses remains approximately constant up to at least 
z^2 at a value a=-1.2±0.1 (consistent with the mod- 
els in Nagami ne et al.|[2005bO . Only at very low masses 
(A4<10 9 Mq at z<l and A4<10 10 M Q at higher red- 
shifts), the SMF seems to become steeper (based on our 
results and those from other surveys) , but this steepening 
has a minor effect on the global stellar mass density. 

5.5. The evolution of the cosmic stellar mass density 



The SMFs were integrated for all masses above the 
completeness level to obtain the observed cosmic co- 
moving stellar mass density. We also integrated the 
Schechter fits to estimate an extrapolated value of the 
cosmic stellar mass density at each redshift interval. In 
Figure O we present these results, comparing them with 
other estimations of the stellar mass density available in 
the literature (see the captions of Figures 0] and [5] for 
references). Note that the observed density values are 
very similar to the extrapolated ones up to z^2, i.e., 
our survey is detecting most of the galaxies that domi- 
nate the global stellar mass density at z<2. We calcu- 
late field-to-field variations of the stellar mass density of 
the order of 30%-40% (depending on the redshift) with 
respect to the average density. For example, the stel- 
lar mass density locked in galaxies with A^>10 11 Mq 
at z>2 is 15%-20% lower in CDF-S than the average 
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of the 3 fiel ds, a slightly lower under- density than that 
observed by Ivan Dokkum et all (|2006l ) comparing three 
^100 arcmin 2 fields (they calculate a 40% difference of 
the CDF-S with the average). 

Figure [5] shows that there is a relatively large increase 
(by a factor of ~1.4) in the stellar mass density of the 
Universe in the last 4 Gyr (from z~0.4 to z=0). This 
large difference could be due to an overestimation of the 
local stellar mass d ensity (suggested by, for example, 
iFontana et all 120041 ) or an underestimation of the den- 
sity at z^0.3 (for example, if low mass objects below our 
detection limit have a non-negligible contribution to the 
stellar mass density at this redshift). However, all the 
estimations of the local density are very similar (diffcr- 
en ces of less than 5% between our value a nd t hose found 
bvlSalucci fe Persid 19991 1 Cole et al.ll200ll. andlBell et al 



2003; even higher values are f ound by [F ukugit a et al. 
19981 iKochanek et alJl200H and lGlazebrook et alj |2003). 



and the sa me occurs for the differen t es timations at 
0.2< z<0.4 (|Brinchmann fc Ellisl 120001 and iBorch et all 
2006). As we discussed in Section 15.41 this significant 
recent evolution of the stellar mass density is mainly 
due to a ^60% increase in the number density of galax- 
ies with 10 9 < M^IO 11 Mq. Assuming an aver- 
age value of the cosmic SFR density of approximately 
0.03 X pyr" 1 at 0.0<z<0 .4 (|Hopkins fc Beacoml 120061 
see also lTresse et al.lfeOOTl and Figure [7]), and a 28% gas 
recycle factor (see Section \TA\ for details), we calculate 
that the stellar mass density of the Universe has grown in 
10 s.o±o.i Mq Mpc -3 from z= o.4 to z=0.0 (in ~4.3 Gyr) 
by just star formation. This is 55±10% of the stellar 
mass density change at z<0.4. Therefore, the remaining 
change in stellar mass density (~10 7 - 9±OA .M0 Mpc~ 3 ) 
must have occurred by either accretion of small satel- 
lite galaxies or by major mergers between gas-depleted 
galaxies (i.e., mergers accompanie d by very lit t le sta r 
formation), as also suggested by iTresse et al.l (|2007h . 
In addition, given that both in the local Universe and 
at 0.2<z<0.4 the SMFs steepen at low stellar masses 
(Ai<10 9 - Mq), the minor merger possibility (accretion 
of A^<10 9 Mq galaxies producing very few or even no 
new stars at all) seems to be favored in detriment of the 
existence of major mergers. 

The evolution in the previous 3-4 Gyr (between z~1.0 
and z~0.4) was slightly slower. About 25% of the lo- 
cal stellar mass density was assembled in that period, 
adding up a total decrease of about 50% in the stellar 
mass density from z=0 to z=l. 

At z~1.0 (8 Gyr ago), the evolution of the stellar mass 
density of the Universe becomes faster (approximately a 
factor of 2), just when the cosmic SFR density re aches 
a maximum (see, e.g.. iPerez-Gonzalez et ai]|2005l ) and 
the galaxies with .M>10 105 Mq dominate the produc- 
tion of stars in the Universe. The rate at which the 
Universe is creating stars stays at approximately a con- 
stant level or decays very slowly from z^l up to at least 
z^ 2 (10 Gyr ago). Between z~l and z~2, the density of 
galaxies with yV(>10 105 Mq decreases significantly (by 
a factor of 3-4). This population of galaxies evolving 
rapidly at l<z<2 (in about 2 G yr) seem to be domina ted 
by early-type objects (see, e.g.. lAbraham et aT]|2007l ). 

Beyond z^2, the errors in the stellar mass density es- 
timates and the differences between the observed and 
extrapolated values become increasingly larger. We find 



that the rate at which stars are being formed remains 
constant or even increases slightly, while the giant galax- 
ies with A4~10 120 Mq are finishing the assembly of 
most of their stellar mass. 

These different steps in the assembly of the cosmic stel- 
lar mass density depicted in Figure [5] are consistent with 
the latest results on the evolution of the observed U V lu- 
minosity density of the Universe (iTresse et al.ll2007D and 
the evolution of the SFR dens ity (jPerez-Gonzalez et al.l 
120051 : iHopkins fc Beacomll2006f) up to z~5. The luminos- 
ity density presents a maximum at around z=1.2 with a 
value approximately 6 times larger than the local UV lu- 
minosity density. At z>1.2, the luminosity and SFR den- 
sity evolution is consistent with a constant. Our results 
are also consistent wit h the hydrodynamical models of 
Na gamine et al.l (|2006f ). which predict that -60% of the 
present stellar mass density was already formed by z=l. 
However, the discrepancy is significant at z>l, where 
these models predict a larger stellar mass density than 
any observation (i.e., they predict a quicker formation 
of the most massive galaxies). The semi-analytic models 
of ICole et all (|2000l) match our results better at z=3 4, 
where they predict a stellar mass density of about 10% 
the present value, but they fail to reproduce the evolu- 
tion at low redshift. 

5.6. Quantifying "downsizing" 

The previous discussion about the evolution of the 
cosmic stellar mass density is clearly consistent with a 
"downsizing" scenario for galaxy formation. We quan- 
tify some properties of this "downsizing" theory in Fig- 
ure [H were we plot the fraction of the total local stellar 
mass density already assembled in galaxies of a given 
stellar mass at each redshift. This Figure shows that the 
most massive systems (.M >10 12 ' Mq, orange widest 
continuous line) formed first (they assembled more than 
80% of their total stellar mass before z=3) and very 
rapidly (about 40% of their mass was assembled in 
1 Gyr between z=4 and z=3). Systems with masses 
lO n - 7 <A4<lO 12 O Al0 assembled their stellar mass more 
slowly: from z— 4 to z— 2.5 (1.5 Gyr), they assembled 
around 50% of their stars, and then evolved more slowly 
to reach the local density at low redshift. Less massive 
galaxies assembled their mass at even a slower speed, 
reaching the local density at very recent epochs. Again, 
this plot shows the rapid recent evolution of the galaxies 
with masses A4~10 10 5 Mq, which assembled ~30% of 
their mass in the last 3 Gyr. 

Our results are c onsistent with the s tellar population 
models assumed by iBrown et al.l (|2007ft for the most lu- 
minous (and probably most massive) red galaxies. Ac- 
cording to that paper, red massive galaxies start form- 
ing at an early epoch, at z=4, followin g an exponen- 
tial S F law with a short r=0.6 Gyr. I Jimenez et al.l 
(2007) find that the most massive early-type galaxies in 
the local Universe formed at z>2.5 and experienced a 
very rapid chemic a l enr ichment, lasting 1-2 Gyr. Also 
van der Wei et all (l2005h find signs of the formation 
of massive (M>2 x 10 11 Mq, according to these au- 
thors) early- type galaxies at z>2, while less massive sys- 
tems present lower formation redshifts (l<z<2). The 
analysis of optical spectra for sphe roidal and bulge- 
dominated galaxies at 0.2<z<1.2 by iTreu et alj (2005) 
also reveals that most of the mass (99%) in systems with 
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Fig. 5. — Evolution of the stellar mass density of the Universe as function of redshift (shown with a logarithmic scale in the quantity l+z 
in the bottom horizontal axis and the corresponding look-back times in the top axis). Our estimations are plotted with black filled stars 
(based on the integration of the stellar mass functions with a Schechter parametrization) and open circles (observed values down to the 
completeness level). Color points and error bars show other estimations found i n the literature. To the refe rences mentioned in the c aption 
of Figure g] we have als o added estimations fromlSahicci fc Persic] !fl999l. SP991 . IBrinchmann k Ellis! pOOOl . BE001 . [Bell et al.1 puM B03), 
Dickin soneTaTI d2003bl . D03), IGlazebrook et al.l 112004 G04), and Rudni ck et al.l [12009 , R06). The inset shows the same evolution of the 
stellar mass density of the Universe, but this time with a linear scale in look-back time in the horizontal axis. 



.M>10 115 Mq formed at z>2, while most recent bursts 
(at z~1.2) can account for 20%-40% of the total stellar 
mass of galaxies with A4<10 110 Mq. IGla zebroo k et all 
(2004) estimates that 38±18% of the stellar mass density 
in galaxies with M>10 108 Mq were already in place 
at z=l, consisten t with our value of ~45%. At z=1.8, 
IGlazebrook et all (|2004h obtain 16±6%, also m agree- 
ment with our own estimation of ~21%. 

If we consider a high value of the fraction of the stel- 
lar mass density already assembled at a given redshift, 
above which the star formation in a galaxy should be 
relatively low, e.g., 70%, Figure [6] shows that galax- 
ies with .M~10 10 - 5 Mq reached that level at z~0.2, 
systems with .M~10 1125 Mq at z~0.4, and galaxies 
with M~10 1175 Mq around z~0.7. These numbers are 
rough ly consistent with the (1+z) 3 5 evolution estimated 
by iBundv et al.l ((2006) for the quenching stellar mass 
(Mq), a mass limit above which the star formation ap- 
pears to be suppressed. 

According to Figure® 50% of the stars in galaxies with 



M>10 110 Mq were already in place at z~0.9. This 
compares well with t he prediction from the models in 
iDe Lucia et al.1 (|2006| ) which establish that half of the 
stars in objects of this mass are assembled into single ob- 
jects at z~0.8. However, these models also predict that 
most of these stars were already formed at z^2.5, but 
were placed in several objects that would coalesce into 
a single object later in a hierarchical way. Our results 
favor a dual scenario where the most massive systems 
with A^>10 12 ' Mq formed most of their stars at z>2.5 
and assembled very rapidly in a way closer to a mono- 
lithic collapse than to a hierarchical coalescence. At the 
same time, less massive systems (10 11 '°>A'1>10 12 ' Mq) 
could have formed their stars later and/or assembled half 
of their mass from several progenitors (where stars were 
already formed at z~2.5) in the time interval between 
z~2.5 and z~l (about 4 Gyr), and most of their mass 
(80%) not before z~0.5. To confirm this scenario, it 
would be necessary to probe the stellar mass function 
at low masses (for objects that would act as building 
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sembled at a given redshift for several mass intervals (wider lines 
referring to more massive systems). Only results for masses above 
our 75% completeness level at each redshift are shown. 



blocks for the galaxies with 7W>10 110 Mo), but the 
scatter of the currently available SMF estimations at low 
masses in this redshift range is too large to obtain ro- 
bust results (maybe due to cosmic variance effects). In- 
deed, our estimations of the cosmic stellar mass density 
at z>3 are affected by the large uncertainties at masses 
below 10 11 M.q (this explains the large difference be- 
tween the observed and extrapolated values of the den- 
sity at z>3). The dual galaxy formation scenario (quasi- 
monolithic and rapid collapse of the most massive galax- 
ies which cease to form stars at a certain epoch, and 
hierarchical collapse for less massive systems) has been 
reproduced by other models where AGNs are supposed 
to q uench the star formation in very mas sive halos (see, 
e.g.. lCroton et al.ll2006t iBower et al.ll2006[ K 



6. 



THE NATURE OF THE IRAC SAMPLE: COMPARISON 
WITH OTHER SURVEYS 



In this Section, we will discuss the main properties of 
the sources in our IRAC sample, comparing them with 
the populations of galaxies detected with different selec- 
tion techniques by other surveys. The results discussed 
in this Section are summarized in Tables [3] and [4] 

Based on the observed photometric data points and 
the SED fit for each galaxy in our sample, we es- 
timated synthetic observed magnitudes in 9 bands 
(FUV, NUV, U n , B, G, K, z, J, and K„) in or- 
der to test which of our galaxies would qualify as Ly- 
man Break Galaxies (LB Gs) at z=1.5-2.5 (LBG-BM 
and LBG-BX galaxi es in ISteidel et aD 120041) . at z~3 
("classical" LBGs. ISteidel et al.1 l2003h. and at z~l 
(GALEX LBGs, IBurgarella et all I2006D . and which 
of our galaxies would qua l ify as Distant Red Gala xies 
(DRGs; iFranx et alJ 120031: Ivan Dokkum et all 120031 ) or 
BzK sources dDaddi et al.ll2004f). Ou r an alysis is similar 
to that used bv lQuadri et al.l ( 20071 ) and IGrazian et all 



<2QQ2b. 

We identifi e d LB Gs following ISteidel et al.l (|2003l ) and 
ISteidel et all (|2004f) . which establish the locus of LBGs 
in a Un — Q vs. Q ~ 1Z color-color diagram, and adopt 
the magnitude cut (i?<25.5) for their survey. We iden- 
tified GALEX LBGs with an analog procedure, but this 
time using a color c riteria based on GAL EX UV pho- 
tometric ba n ds (se e IBurgarella et ail 120061 ). Following 
IFranx et al.l ([2003), we defined DRGs as the galaxies 
presenting a color J — K s >1.37 [corresponding to (J — 
-Kg)vega>2-3]. Finally, we identified star-forming BzK 
galaxies (BzK-SF) and passively evolving BzK galaxies 
BzK-PF) using Equations (2) and (3) in Dad di et al.l 
2004). In our IRAC selected sample, we identified 6,656 
sources as LBGs with i?<25.5 (summing up all types), 
763 sources as DRGs with K<22.9, and 2,426 as BzK 
sources with K<22.9 and z>1.4 (summing up the two 
types). 

The average surface density of LBGs (including all 
sub- types) with i?<25.5 detected by our IRAC survey 
is 10.0 LBGs arcmin -2 . We detect 0.7 LBGs arcmin -2 
with the GALEX bands and iVW<25.0, a higher 
density than the one given by IBurgarella et all (2006, 
0.3 arcmin -2 ), but closer to the density given in 
IBurgarella et~all (|2007l . 1.0 arcmin" 2 ). We find 4.6 LBG- 
BMs arcmin -2 (5.3 LBG-BMs arcmin -2 without any op- 
tical magnitude cut), 3.1 LBG-BXs arcmin - (3.6 LBG- 
BXs arcmin -2 without any optical magnitude cut), 
and 1.6 "classical" LBGs arcmin -2 (2.0 "classical" 
LBGs arcmin - at any /-^-magnitude) , very similar val- 
ues to those found in ISteidel et ait (|2004l 5.3 LBG- 
BMs arcmin -2 and 3.6 LBG-BXs arcmin -2 ) and 
Steidel et all d2003l 



1.7 "classical" 



arcmin 
LBGs arcmin - ^). 



Adelberger et al.l (|2004) and IGrazian efaD (2007) also 



find very similar surface densities for the 3 types of LBGs 
at z>l. The median magnitudes for the LBG sub-sample 
are i?=24.6 and iC=23.1, a very faint NIR magnitude 
only reachable by the deepest ground-based or IRAC sur- 
veys. The average photometric redshifts for the LBGs 
in our sample are consistent with the literature, as dis- 
cussed in Appendix [HI jointly with the results for the 
other populations of high redshift galaxies. 

The average surface density of DRGs in our IRAC 
survey is 1.8 DRGs arcmin -2 (1.1 DRGs arcmin -2 for 
sources with i\~<22.9), a value in between the densi- 
ties quoted bylFranx et"ail (|2003l . 3.0 DRGs arcmin -2 ), 
iForster Schreiber et all (|2004L 1.0-1.6 DRGs arcmin -2 ), 
and iPapovich et al.l (|2006t 0.8 DRGs arcmin -2 ). The 
median magnitudes for the DRG sub-sample are i?=25.7 
(a very faint optical magnitude beyond the reach of most 
UV/optical surveys) and K—22.6. 

The average surface density of BzK galaxies at 
z>1.4 down to K=21.9 is 1.3 BzKs arcmin -2 , di- 
vided into 0.2 BzK-PF arcmin -2 and 1.1 BzK- 
SF arcm in -2 . This i s con sistent with the densities 
given in iDaddi et all (|2004 0.22 BzK-PF arcmin -2 
and 0.91 BzK-SF arcmin -2 ) for the same brightness 
limit. At fainter magnitudes, K<22.9, we identify 
0.4 BzK-PF arcmin -2 and 3.3 BzK-SF a rc min -2 
close to the values found by iReddv et al.l (|2006bl 
0.24 BzK-PF arcmin -2 and 3.1 BzK-SF — ; - -2 



arcmin 



) 

and IGrazian et al.l (|2007l 0.65 BzK-PF arcmin 2 and 
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3.2 BzK-SF arcmin -2 ) for the same magnitude cut. At 
even fainter if -band magnitudes, the source density of 
galaxies identified as BzK continues rising (especially 
the SF sub-type) as redshift interlopers become more nu- 
merous (up to 40%). 

These figures (densities and average redshifts) demon- 
strate that our IRAC survey constitute an almost com- 
plete census of the previously detected galaxies at 
1.5<z<4, including most of the LBGs, DRGs, and BzK 
sources, the most important populations of galaxies se- 
lected at z>l. Still, some of the IRAC sources are not 
recovered by any of these selection criteria (even when no 
magnitude cut is performed for LBGs, DRGs, or BzK 
sources). The numbers of these galaxies recovered only 
by the deep IRAC observations are given in Table [3l 

The LBG population accounts for a negligible fraction 
(less than 10%) of the entire IRAC sample at 0.4<z<1.0. 
At z<0.4, ~30% of IRAC galaxies are classified as LBGs 
(40% at 0.0<z<0.2 and 20% at 0.2<z<0.4), most of them 
within the LBG-BX sub-type, which has a significant 
fract ion of z<l interlope rs at bright apparent magnitudes 
(see ISteidel et all 12004) . LBGs selected with GALEX 
bands are also a minor fraction (around 5%) of the to- 
tal number of IRAC galaxies at 0.8<z<1.3. However, 
at z>l, other LBG sub-type start to be very numerous 
and even dominate the IRAC galaxy counts: ^35% of all 
the sources in our IRAC survey at 1.0<z<1.3 are LBGs 
(80% of them LBG-BMs), 50%-60% at 1.3<z<3.0 (with 
similar contributions from the different sub-types), 65% 
at 3.0<z<3.5 (all of them "classical" LBGs), and 50% at 
3.5<z<4.0 (all of them "classical" LBGs). These frac- 
tions are slightly higher for LBGs not limited by any 
i?-band magnitude. 

The median stellar masses of LBGs range from 
10 9 - 6 M Q to 10 10 - 2 M Q at l<z<2.5. These values are 
0.1-0.2dex lower than the median stellar masses for the 
global population of IRAC sources at each redshift in- 
terval. For this reason, although their numbers are rel- 
atively large (even dominate the number counts), LBGs 
have a less important contribution to the global stel- 
lar mass density. Indeed, at 1.0<z<1.6, they harbor 
less than 25% of the total stellar mass density 6 . At 
1.6<z<4.0, they account for 35%-45% of the total stel- 
lar mass density (ro ughly consistent with th e estimations 
from the models in lNagamine et aT1l2005bl ). These per- 
centages increase by 5%-15% if we consider all LBGs 
without any i?-band c ut, then making our estimations 
consistent with those in lGrazian et al.l (|2007l . where they 
do not apply any magnitude cut). 

Around 10% of LBGs are detected by MIPS at 24 ,um 
above the 75% completeness flux level (20% with any 
flux), especially at 1.6<z<2.5, where MIPS is more effi- 
cient detecting sources due to the pass of the 7.7 /xm PAH 
feature through the filter. At z>2.5, the fra ction of MIPS 
detect ions is about 10%, consistent with iHuang et al.l 
(|2005t ). MIPS detections are more common for the high- 
est mass galaxies: the median stellar mass for 24 /an 
detected LBGs (10 109 M Q at z-2.5 and 10 111 M Q at 
z~3.0 for sources with F(24)=80 /xJy) is ~0.8dex larger 

6 This percentage has been calculated by adding the total stellar 
masses of LBGs in that redshift interval, and dividing it by the total 
stellar masses of all galaxies. This must be analogous to dividing 
the stellar mass densities of both galaxy populations for a fixed 
volume (that enclosed at the given redshift interval). 



than the median for all LBGs. LBGs with MIPS detec- 
tions account for 10%-20% of the total stellar mass at 
z>1.5. 

In contrast with the previous figures for LBGs, DRGs 
are less numerous but more massive. DRGs only ac- 
count for 15% of the sources at 2.0<z<2.5, and ~30% 
at 2.5<z<4.0. However, their median stellar masses are 
larger than those of LBGs at each redshift, and even 
larger than the median for the entire population of IRAC 
sources. For l<z<2, their median masses are 0.3-0.5dex 
larger than those for the entire IRAC population, and 
at z>2 they remain ^0.6dex larger (with a median of 
10 110 A4q). This translates to DRGs accounting for 
70% of the total stellar mass density at z>2.5, ~35% 
at 2.0<z<2.5, and less than 10% below z =2. These fig- 
ures a re very similar to those found by Rudnick et al. 
(2006), who find that DRGs contribute 30% and 64% 
to the stellar mass density at z~2 and z—2.8, respec- 
tiv ely. They are a ls o con si stent with the resu l ts obt ained 
bv IGrazian et al.l <l2007T) . iMarchesini et al.l (|2007f ) and 
Ivan Dokkum et all (|2006f) . Note that most of the stellar 
mass density of the Universe at z>2.5 would not be de- 
tected by optical surveys re aching depths br i ghter than 
i?~25.5. C o nsiste ntly with iPapovich et al.1 (|2006f ) and 
IWebb et al] (|2006l) . we find that within the DRG popu- 
lation, about 40% are detected by MIPS at 24 (up to 
50% at 2.0<z<3.5), and these objects have median stellar 
masses 0.1-0.3dex larger than the median for all DRGs. 
DRGs with MIPS detections account for more than 40% 
of the total stellar mass at z>2.5 (20% at 2.0<z<2.5, and 
less than 5% at z<2.0). 

The BzK criterium is very effective in detecting mas- 
sive galaxies at z>1.5, even more than the J — K selec- 
tion of DRGs. Up to 75%-95% of the IRAC sources at 
1.3<z<2.5 are recovered by the BzK selection, 80% at 
2.5<z<3.0, 55% at 3.0<z<3.5, and 30% beyond z=3.5. If 
we only consider BzK galaxies with K<22.9, these per- 
centages decrease by a factor of ~2. Most of the BzK 
galaxies are classified as star- forming (typically 90%). 
Median masses for BzK galaxies range from 10 10 2 at 
z=1.5 to 10 10 - 6 at z=3.0 and 10 109 at z=4.0, 0.1-0.4 dex 
larger than median stellar masses for the whole IRAC 
sample. This translates into BzK galaxies tracing a large 
fraction of the stellar mass density at z>1.5: more than 
55% and up to 97% at z>1.3. Again, if we only consider 
BzK galaxies with K<22.9, these percentages decrease 
by 15%-20%. These fractions are comparable to the 94% 
contribution of BzK sour ces to the total stell ar mass 
density at z— 1.8 found bv IGrazian et all (|2007[ ). Typi- 
cally, 30% or more BzK galaxies are detected by MIPS at 
24 (Uirx, with a predilection for the passively evolving sub- 
type at z>2 (~60% and ~30% of BzK-PE and BzK-SF 
galaxies are detected by MIPS). This means that pas- 
sively evolving BzK galaxies may still harbor significant 
star formation or obscured AGNs. 

Very few galaxies are identified as LBGs and DRGs si- 
multaneously in our IRAC survey: just 5% of all galaxies 
at 2.0<z<3.0, ~8% at z>3, and less than 1% elsewhere. 
However, this does not mean that the 2 selection crite- 
ria are completely orthogonal. Indeed, about 20% of the 
DRGs at z=2-3 and 30% of the DRGs at z>3 qualify as 
LBGs, and 15% of LBGs at z>2.5 are DRGs. Most of the 
LBGs that also qualify as DRGs lie in the "classical" sub- 
type (more than 95% of them), which makes our results 
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also co nsistent with the fractions found in lGrazian et all 
(|2007ft . who only discussed BM-BX objects. Note that 
if we only consider the DRGs brighter than K=21.9, the 
fraction of LBGs that are also DRGs drops below the 5% 
level, in good agreement with the 10% up per limit predic- 
tion fro m the hydrodynamic models of INagamin e et ail 
(|2005aD . 

The BzK and DRG selection criteria present a large 
overlap. Around 20%-30% of all IRAC galaxies at z>2 
are recovered by both selection techniques, especially by 
the BzK-SF criterium. Indeed, more than 95% of all 
DRGs at 1.5<z<3.0 are BzK galaxies, most of them 
(^90%) wi thin the star-forming BzK sub-type (in agree- 
ment with lGrazian et aT1l2007f ). DRGs are only a minor 
contributor to the BzK population at z<3, where less 
than 35% of BzK sources are DRGs, but this percentage 
rises to 70% at z>3.5. It is also interesting to mention 
that 50% of BzK-PE galaxies at 2.0<z<2.5 and all the 
BzK-PE galaxies at z>2.5 are DRGs. 

None of the LBGs lie in the BzK-PE type (as also 
noted bv lGrazian et al.ll2007T) . but the BzK-SF type also 
has a large overlap with the LBG population: more than 
95% of BM-BX galaxies up to z~2 are recovered i n the 
BzK diagram (consistent with iReddv et al.ll2005[ ). and 
less than 40% in the case of "classical" LBGs at z>3. 

7. LINKING STELLAR MASSES AND STAR FORMATION 
RATES UP TO Z=4 

7.1. The evolution of the cosmic star formation rate 

density 

The time derivative of the stellar mass density function 
plotted in Figure [5] can be used to obtain the evolution 
of the SFR densit y of the Universe , i.e., the well-know n 
Lilly-Madau plot ([Lilly et al.lll996t iMadau et al.lll996f ). 
For each pair of stellar mass density points in Figure [SJ 
we have estimated the SFR density (averaged through 
the time interval enclosed by the points) necessary to 
produce the stellar mass density difference between the 
corresponding redshifts. This SFR density must be cor- 
rected upwards by some amount to account for the mass 
loss due to stellar winds and superno ya ejecta. For a 
Salpeter IMF, the correction is 28% (ICoIe et all 120011 : 
iDickinson et~aill2003bt iHopkins k. Beacomll2006i r These 
SFR density estimations are plotted in Figure [7] with 
filled black stars (open circles show the average SFR den- 
sities derived from the observed stellar mass density val- 
ues), and compared with other cosmic SFR density es- 
timations (based on direct SFR measurements) found in 
the literature. Surprisingly, our estimations of the cosmic 
SFR density are systematically smaller than the previ- 
ously published results (on average, a factor of ~1.7 at 
z<2 and a factor of 4.5 at higher redshifts, compared to 
individual es timations). Thi s discr e pancy has also been 
remar ked b y iRudnick et alj (120061 ). IHopkins fc Beacoml 
(|2006f ). and iBorch et all (|2006l )~~who integrate the time 
evolution of the cosmic SFR density to obtain the evolu- 
ti on of the stella r mass density. 

IRudnick et all (|2006h find a factor of 1.6-2.5 offset at 
z<2 between the measured stellar mass densities and the 
values derived from th e integration of t he SFR density 
evolution published bv lCole et~aT1 (|2001h . At z>2, they 
find very good agreement. Ho wever, the fit of t he SFR 
density evolution published bv lCole et alj (|2001l ) gives a 
factor of 2-3 lower SFR densities than the latest estima- 
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Fig. 7. — Evolutio n of the co-movin g SFR density of the Universe 
(Lilly-Madau plot, ILillv et alj 119961; IMadau et alj|l996ft . Filled 
stars and thin error bars show the SFR density estimations based 
on the time derivative of the stellar mass density evolution shown in 
Figure [5] Open circles show the derivative for the observed values 
of the stellar mass density. The colored points (shown with error 
bars) are extracted from different sources in the literature (using 
different S FR trace r s), co mpil ed and normalized to the s ame cos- 
mology by Hopkmi 12004 ) and IHopkins fc Beacoml H2009) . To this 
compilation, we have also added th e SFR density estimations at 
z~2 and z~3 found in Rcddy et al. (2007). Red symbols are es- 
timations based on Hydrogen emission-lines, and green points on 
[0//JA3737 SFR estimations. UV-based data points are plotted in 
blue. Cyan estimations are based on mid-infrared data. Magenta 
points are based on sub-mm and radio observations. The yellow 
point is based on X-ray data. Thick black error bars show weighted 
averages and standard deviations of the literature data points for 
the 12 redshift intervals considered for the stellar mass functions in 
this paper. The black line shows the evoluti on of the cosmic SFR 
density as parametrized in Cole et al. (2001). 



ti ons at z>2 (see Figu r e E). 

Hopk ins & Beacoml ((2006) argue that the difference 
can be related to a limitation in our understanding of 
the IMF, given that the direct SFR estimations are sen- 
sitive to the high mass end of the IMF, while the stellar 
mass estimations are sen sitive to the low mass end. In- 
deed, IBorch et al.l (|2006l ) find good agreement between 
the SFR and the stellar mass density evolution at z<l 
by choosing a lChabrierl (|2003[) IMF, w hich gives stellar 
masses similar to a lKroupa et al.l (|1993l ) IMF (~1.7 times 
lower than our masses, based on a Salpeter IMF), but 
with lower SFRs (by a factor of 3). In our case, an offset 
of 3/1.7~1.8 would also make consistent the results of the 
evolution of the SFR density and s tellar mas s densi ty up 
to z~2. However, at z>2, the same lChabrien (|2003D IMF 
fails to match the SFR and stellar mass densities: using 
the same Chabrier IMF, a good fit to the stellar mass 
density evolution gives SFR densities lower than the lat- 
est observations by up to a factor of 2-3. A top-heavy 
IMF (compared to a Chabrier IMF) at high redshifts (i.e., 
an evolution of the IMF) could make the SFR and stellar 
mass density evolutions at z>2 match. A top-heavy IMF 
has also been proposed by the galaxy evolution mod- 
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els in iBaugh et all (|2005f ) and lLacev et all (|2007h to ex- 
plain the number counts of galaxies in the I R and sub- 
mm spectral ranges (see also , among others. |Elmegreerjj 
200 i iNaeashima et a.l.ll2005t ISiigi2005t iLe Delliou et all 
200a iKlessen et al.ll2007f ). Theoretical calculations of 
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formation of stars also predict a top-heavy IMF for star- 
burst s (e.g., iPadoan et~aTl Il997t iLarsonl [l998t iKroupal 
2007). Finally, some observational evidence also sug- 
gests a top-heavy IMF for certain stellar populations 

Ivan Dokkum &: van der Marell 



( Franceschini et al 



2007; M aness et al 



2001 



2007) 



The discrepancy in Figure [7] could also be solved if 
the SFRs estimated (with different star formation trac- 
ers, mainly the UV, IR, and sub-mm emission) for the 
massive galaxies at z>2 were overestimated due, for ex- 
ample, to the presenc e of strong obscure d AGN in most 
of these sources (see iDaddi et al.l 120071 and references 
therein), which would imply that a significant fraction of 
their UV or IR emission arise from the AGN, i.e., it is 
not linked to star formation. 

Finally, the reader should also note the very high 
SFR density derived for 0.0<z<0.4 from the stellar mass 
density derivative, directly related to the significant in- 
crease in the stellar mass density of the Universe in 
this time period, as discussed in previous Sections. Di- 
rect SFR density measurements at z<0.4 are a factor of 
2-4 smaller, which indicates that the evolution of the 
stellar mass function at low redshift is not only gov- 
erned by star formation (but also by mergers) and/or 
we may be underestimating the stellar mass density at 
z^0.4 (if there is a numerous population of low mass 
galaxies below our detection limit, which can merge to- 
gether in the last 4 Gyr to increase the density of galaxies 
with mass A^IO 9-11 M.q) or overestimating the local 
value. In Section 15.51 we estimated a change of the stel- 
lar mass density of 10 7 - 9±01 M G Mpc~ 3 from z=0.4 to 
z=0.0 due to dry mergers. Since these mergers are af- 
fecting galaxies with M — IO 9 ^ 11 A4q and the average 
number density of these systems is 10~ 2,1 Mpc~ 3 log(.M) 
at 0.2<z<0.4, we calculate a mass accretion rate of 
AM/M=0.12±0.03 Gyr" 1 from z=0.4 to z=0.0. This 
value is in good agreement with the accretion rate 
of z=0.1 gal a xies in the red sequence, estimated by 
van Dokkuml (j2005h in 0.09±0.04 Gyr" 1 (our value is 
just 30% higher, but consistent within errors), who also 
calculate that the median mass ratio of the mergers in 
nearby early- type galaxies is 1:4. 



7.2. Specific star formation rates 

Figure [5] shows the evolution of the specific SFR (SFR 
per stellar mass unit) of the Universe (gray shaded area), 
calculated by dividing the average values of the cosmic 
SFR density (given in Figure |7|) by our stellar mass den- 
sity estimates (given in Figure [5]) . There is a continuous 
increase of the specific SFR of the Universe as we move to 
higher redshifts. If we consider the evolution of the spe- 
cific SFR for different stellar mass intervals (color lines in 
Figure |8|) , we clearly see that the most massive galaxies 
(A4>10 117 M-q) presented very large specific SFRs at 
high redshift. These galaxies exhibit values of the SFR 
that are so large that they could double their stellar mass 



o 



Li_ 



o 



10 8 '°< M/M < 1CT 
10 9 -°< M/M <U)" 
10 m0 < M./M B <10 U 

10 1, < M/Jl s <10 U J 
M./M @ >W r 




Fig. 8. — Evolution of the cosmic specific SFR. The black con- 
tinuous line shows the evolution obtained from the the stellar mass 
and SFR density estimates plotted in Figures [5] and [7] The gray 
shaded area depicts the typical uncertainties in the calculation of 
the cosmic specific SFR. Color lines join the median values of the 
distribution of specific SFRs of our IRAC sample for several mass 
intervals (excluding all X-ray detected sources) , while vertical error 
bars show the quartiles of that distribution. For clarity, we have 
only depicted the quartiles for non-consecutive mass intervals. We 
only show the median and quartiles for redshift bins where we de- 
tect more than 10 galaxies within a given stellar mass interval. 
Dashed lines mark the redshift ranges where our sample is less 
than 75% complete for the given mass interval. 



in just 0.1 Gyr (see the scale on the left axis) at z=3-4 7 . 
As we move to lower redshifts, their specific SFRs de- 
crease considerably, by a factor of 10 from z^4 to z~2.5 
(in less than 1.5 Gyr), and by a factor of ~100 from z~4 
to z^l.5 (in 3 Gyr), i n agreement with the results from 
iPapovich et al.l (|2006f) . For lower stellar masses, the evo- 
lution is less pronounced. For example, galaxies with 
10 10 <A / J<10 110 Mq, whose evolution is very similar 
to the cosmic average, present a decrease in the specific 
SFR of a factor of 10 from z~2.5 to z-0.5 (in 6 Gyr). 
The evolution at z<l is very simil ar for all the stella r 
mass intervals, as already noted by I Zheng et al.l (|2007|) . 
with a change of about a factor of 10 in this period of 
~8 Gyr. In this period, there is a significant luminosity 
(and maybe density) evolution of luminous IR galaxies 

7 The individual SFRs and median specific SFRs for the entire 
sample derived from the extrapol ated rest-frame 24 )im lumin osi- 
ties and the calibration given by Alonso- Herrero etall (2006, see 
also ICalzctti et al. 2003) are a factor of 2-3 smaller at z>2 than 
those derived directly from extrapola ted estimations of the total 
IR luminosity. Papovich et al. (2007) also find (based on 70 /im 
observations) that SFRs of z~2 galaxies obtained from total IR lu- 
minosities (extrapolated from fits of dust emission models to single 
24 detections) are overestimated by factors of a few. Using 
the SFR estimations based on rest-frame 24 fim luminosities, the 
median specific SFRs of the most massive galaxies are very similar 
to the Universe average (the shaded area in F igure [SJ at z=3— 4, 
and our global results presented in Section l7.2l are not affected. 
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(|Charv & Elbazll200ll : IPerez- Gonzalez et al.ll2005f) . This 
evolution is also detected in our stellar mass analysis, 
given that the fraction of the total stellar mass density 
locked in galaxies emitting strongly in the thermal IR 
(and being detected by MIPS) increase from about 10% 
at z=0 to ~50% at z=0. 7-1.0, where LIRGs dominate 
the cosmic SFR density. The fraction of the total stellar 
mass density locked in MIPS sources remains approxi- 
mately constant from z~l up to z=4, where LIRGs and 
even ULIRGs have a significant contribution to the total 
SFR density, as the comparison w ith other SFR tracers 
show (Perez- Gonzalez et alll2005h . 

The evolution of the cosmic specific SFR is also con- 
sistent with the "downsizing" picture described in Sec- 
tion 15 -61 where the most massive galaxies formed most 
of their mass at z>3 in very intense and rapid episodes 
of star formation, presenting high specific SFRs which 
would double the stellar mass of these systems in time 
scales shorter than 1 Gyr. Less massive systems assem- 
bled more slowly, presenting specific SFRs which would 
double their mass in time scales comparable to the look- 
back time of the Universe at each redshift. 

8. CONCLUSIONS 

We characterize the mass assembly of galaxies in the 
last 12 Gyr (~90% of the Hubble time) by analyz- 
ing the stellar mass functions and densities estimated 
from a sample of ^28,000 sources selected in the rest- 
frame near-infrared. The sample has been built from 
Spitzer/IRAC (the selection being made at 3.6 /xm and 
4.5 /itm) observations of 3 fields: the Hubble Deep Field 
North, the Chandra Deep Field South, and the Lockman 
Hole. The total surveyed area is 664 arcmin 2 . This IRAC 
sample is 75% complete down to 1.6 /xJy ([3.6]=23.4), 
which translates to an approximate stellar mass com- 
pleteness level (for passively evolving galaxies formed at 
z=oc) of at least 10 9 M® up to z=0.5, 10 10 M® up to 
z=1.0, and 10 11 M® up to z=4.0. In order to analyze 
the effects on our results of low mass galaxies faint at 
rest-frame near-infrared wavelengths, we complement the 
IRAC survey with an optically (I<25.5) selected sample 
of a similar size. 

We estimate photometric redshifts, stellar masses, and 
star formation rates for all galaxies using a set of tem- 
plates built by modelling the stellar population and dust 
emission of galaxies with known spectroscopic redshifts. 
The quality of our photometric redshifts is very good for 
more than 85% of the sample, and good for nearly 95%, 
according to the comparison with spectroscopic redshift 
at z<1.5. At z>1.5, where spectroscopic redshifts are 
scarce, we test our photometric redshifts by comparing 
the average values and standard deviations of the red- 
shift distributions of several galaxy populations: Lyman 
Break Galaxies (LBGs), Distant Red Galaxies (DRGs), 
and BzK galaxies. We find very good agreement with 
spectroscopic surveys of these sources, and other photo- 
metric redshift analysis. We also analyze the goodness of 
the stellar mass and star formation rate estimates, find- 
ing that they are accurate within a factor of 2-3 (see 
Appendix [B]) . 

Our estimation of the local stellar mass function is 
in good agre ement with prev i ous estimations based on 
2MASS data dCole et al.ll200lL iBell et al.ll200l . We find 
a slope of a~-1.2 at low stellar masses (similar values 



are also found for stellar mass functions at all redshifts 
up to z=4), and a pronounced steepening of the stellar 
mass function at A / l<10 9 M®. Approximately 1 out of 4 
local galaxies are actively forming stars. Around 10-15% 
of the global stellar mass density in the local Universe is 
found in active star-forming galaxies (in agreement with 
iPerez-Gonzalez et~a l. 2003a), and this percentage rises to 
~50% at z~l, remaining approximately constant beyond 
that redshift. 

Our results indicate that the most massive systems 
(M>10 12 M®) assembled the bulk of their stellar 
mass in a very rapid collapse (half of their stellar mass 
in less than 1 Gyr) at early epochs (z>3 or 11 Gyr 
ago), close to what can be regarded as a monolithic col- 
lapse. The formation was characterized by large spe- 
cific SFRs with doubling times of about 0.1 Gyr. Galax- 
ies with 10 115 <A^<10 12 M® formed more slowly, as- 
sembling half of their stellar mass before z~1.5 (more 
than 9 Gyr ago) and more than 90% of their stel- 
lar mass beyond z~0.6. Less massive systems (with 
10 9 <.M<10 110 Mq) formed at even a slower speed: 
half of their stellar mass was assembled beyond z~l 
(more than 7 Gyr ago), and they experienced a signif- 
icant increase in their stellar mass (20%-40%) recently 
(at z<0.4 or in the last 4 Gyr), probably by dry accre- 
tion of small satellite galaxies (with an accretion rate 
AM/M=0.12±0.03 Gyr" 1 ). The specific SFRs of these 
galaxies evolved (closely to the cosmic average) from 
10 Gyr -1 at z^4 to less than 0.1 Gyr -1 in the local 
Universe. 

We find that approximately half of the local stel- 
lar mass density was already formed at z~l (8 Gyr 
ago), which translates to an average assembling rate of 
0.048 A^ Q yr~ 1 Mpc -3 (taking into account a 28% recy- 
cle factor, i.e., the fraction of the total stellar mass re- 
injected in the interstellar medium in the form of stellar 
winds and supernova ejecta). At least another 40% of 
the local stellar mass density assembled from z=l to z=4 
(in 4 Gyr) at an average rate of 0.074 0yr _1 Mpc -3 . 
We find that the cosmic SFR densities estimated by 
differentiating the evolution of the cosmic stellar mass 
density do not match the observ ations based o n direc t 
SFR tracers, as als o noticed by iRudnick et al.l (|2006h . 
iHopkins fc Beacoml (|2006h . and lBorch etafl |2006j). The 
mismatch up to z~2 (a fac tor of ~l- 7 ) coul d be explained 
by c hanging the IM F to a lChabrierl (20031 ) IMF (instead 
of a lSalpetenll955l IMF, the default used in this paper). 
At z>2, the discrepancy is larger (a factor of 4-5), and 
can only be solved if the IMF is top-heavy (i.e., a differ- 
ent IMF at high redshift) and/or if the SFRs of the most 
massive galaxies at z=3-4 (calculated with different star 
formation tracers, mainly the UV, IR, and sub-mm emis- 
sion) are overestimated due, for example, to the presence 
of strong obscured AGN in most of these sources (which 
would imply that a significant fraction of their UV or IR 
emission arises from the AGN, i.e., it is not linked to star 
formation) . 

We confirm th at galaxy formatio n follows a "down- 
sizing" scenario (jCowie et al.l 11996ft . Our results are 
broadly consistent with previous obser yational works 
that c onfirm this fo rmati o n theory (Heavens et al.l 
20041: Uuneau et all 120051 : iGfazebrook et all 12004 



Perez-Gonzalez et al.l 120051 : iFontana et al.~l2006l). and 
with models of galaxy formation (e.g.. iNagamine et al.l 
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120041: ICroton et all 120061: iBower et alj H)06). At low 
rcdshift (z<l), there is also an "upsizing" effect, when 
intermediate mass galaxies (yV(=10 9_11 M.q) increase 
their density by accretion or coalescence of previously 
formed smaller galaxies (producing little star formation) . 

We have also analyzed the nature of the galaxies in our 
sample, comparing them with the populations of sources 
detected with different selection techniques by other sur- 
veys. Based on the measured number densities and red- 
shifts, we conclude that our survey constitutes an almost 
complete census of the different populations of previously 
known galaxies at high redshifts, including most of the 
LBGs at l<z<3.5, most of the DRGs at z>2, and most of 
the BzK galaxies at z>1.4. LBGs dominate the number 
counts of IRAC galaxies at high redshift, being about 
a factor of 2-3 more numerous than DRGs and BzK 
galaxies, but most of the stellar mass density (more than 
50% and up to 97%) at z>2.5 resides in the latter, while 
LBGs account for less than 50% of the total stellar mass 
density. 
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TABLE 1 

Stellar mass functions for. the global and star-forming population of 

galaxies. 
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NOTE. — a In units of A1q. b In units of Mpc 3 (log AI) 
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TABLE 1 

Stellar mass functions for. the global and star-forming population of 

galaxies. 
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TABLE 1 

Stellar, mass functions for. the global and star-forming population of 

galaxies. 
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TABLE 1 

Stellar mass functions for. the global and star-forming population of 

galaxies. 



Redshift range log(M) a logOiRAc) 6 log(</>i_band) 6 log(0 M lPs) 6 
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NOTE. — a In units of M@. b In units of Mpc" 3 (log \ M)~ 



TABLE 2 

Results of the Schechter (1976) fits (including stellar mass densities) to the global and star-forming stellar mass functions. 



GLOBAL STAR-FORMING 



K cdshift 


range 




a 




log(M*) a 




l°g(p*) c 


log(pf s ) c 




«SF 




iog(M* SF r 


log(^ P ) 6 


log(pf F ) c 


log(pf s ) c 


0.0 < z < 


0.2 


-1 


.18±0. 


12 


11.16±0.25 


-2.47±0.22 
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8.75 


-1 
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Note. — 



a In units of Mq. b In units of Mpc 3 (log M) 1 . c In units of A4 Mpc 3 . 
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TABLE 3 

The IRAC sample: Comparison with other surveys. 



Number of sourccs/MIPS detections Number of sources 



IRAC" LBG h DRG a BzK a LBG 6 DRG a BzK a JUST 

tj 

OK 



Redshift 


GALEX 


BM 


BX 


"classic" 




PE 


SF 


DRG a 


BzK a 


LBG 6 


IRAC a ' c 


(0.0,0.2 


1473/242 


0/0 


81/11 


510/94 


2/0 


9/4 


0/0 


468/62 


3 


9 


243 


655 


(0.2,0.4 


1745/375 


0/0 


121/10 


183/30 


5/0 


8/1 


0/0 


318/54 


1 


8 
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1252 


(0.4,0.6 
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0/0 
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29 
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(0.6,0.8 


3953/979 


0/0 


55/2 


6/0 


1/0 


38/9 


2/1 


403/55 


1 
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(0.8, 1.0 
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12/2 


0/0 


78/28 


2/0 
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7 


46 
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3314 


(1.0, 1.3 
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189/66 
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117/25 


1/1 


118/48 


21/5 


1060/130 


13 


87 
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(1.3, 1.6 
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0/0 
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5/2 


61/26 
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2 


58 
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315 


(1.6,2.0 


1640/532 


0/0 


414/105 
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24/8 


55/18 


104/33 


1417/468 


8 


55 


821 


90 


(2.0,2.5 


1404/439 


0/0 


144/39 
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95/29 


231/109 


49/17 


1274/406 


43 


230 
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52 


(2.5,3.0 
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0/0 


4/0 


135/31 


278/82 
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23/15 


677/197 


59 


234 


322 


68 


(3.0,3.5 


558/162 


0/0 


0/0 


1/1 


365/85 


171/93 


15/10 


294/104 


49 
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52 


(3.5,4.0 


529/95 


0/0 


0/0 


4/0 


276/36 


165/59 


33/20 


115/23 


40 


102 


65 


117 



Note. — a Any magnitude. b Magnitude limited to i?<25.5. c IRAC sources not recovered by any other selection criteria (i.e., they are not LBGs, DRGs, or 
BzK galaxies). 



TABLE 4 

The IRAC sample: stellar mass statistics and contribution to the stellar mass density. 



Stellar masses' 1 and percentage of total stellar mass density 
Redshift ALL LBG DRG BzK 



IRAC Any magnitude _R<25.5 Any magnitude i^<22.9 Any magnitude i^<22.9 



(0.0,0.2] 
(0.2,0.4] 
(0.4,0.6] 
(0.6,0.8] 
(0.8, 1.0] 
(1.0, 1.3] 
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(2.0,2.5] 
(2.5,3.0] 
(3.0,3.5] 
(3.5,4.0] 
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10.5 
10.5 
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72% 
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87% 


31% 
53% 
62% 
50% 


97% 
87% 
81% 
56% 


83% 
69% 
64% 
39% 



Note. — 



a Logarithms of the median and quartiles of the distribution of stellar masses in units of [A^©]. 
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APPENDIX 

A. THE MERGED PHOTOMETRIC CATALOG 

This Appendix describes how we selected and measured mult i- wavelength photometry for the galaxies included in 
the IRAC and /-band selected samples. First, we characterize the reduction, detection and photometry procedures 
in the Spitzer images. Then, we outline how we merged this photometry with the fluxes estimated in ground-based 
optical and NIR images. Special details are given for the spectroscopy compiled for our sources. We also discuss the 
methods used to remove stars from our catalogs. Finally, we discuss the presence of Active Galactic Nuclei (AGN) in 
our samples. 

A.l. IRAC and MIPS detection and photometry 

We compiled all the IRAC and MIPS data available in the HDF-N, the CDF-S, and the LHF, including the GTO 
data in the 3 fields, the GOODS data in HDF-N and CDF-S, and the data around the GOODS footprint in the CDF-S 
taken as a Spitzer Legacy Survey (PI: van Dokkum) . All the reduced data (Basic Calibrat ed Data products de livered 
by the Spitzer Science Center) were mosaicked together using the procedure developed by iHuang et al.l (|2004f) . This 
procedure includes pointing refinement, distortion correction, drizzling to a scale half of the original (approximately 

0.6 arcsec/pixel), and correction of detector artifacts (more noticeably, mux-ble eding). 

Detection of sources in the IRAC images was performed with SE XTRACTOR (Bertin fc Arnoutsi n"996). Given that 



the FWHM of the IRAC Point Spread Function (PSF) is 1.8-2.0" (jFazio et al.ll2004af) . and that the PSFs are very 



sharp and stable, almost all sources are point-like in the 4 channels, and objects can be resolved for separations of 
the order of The crowdedness of the our very deep images in the two bluer IRAC bands is very high, mostly 

at 3.6 /im and especially near bright stars, making the deblending of sources hard for automatic procedures such as 
that used by SEXTRACTOR. To alleviate this problem, we detected sources at 3.6 ftm and also (separately) at 4.5 /urn, 
where the depth is slightly lower and crowdedness is less severe. The two catalogs built in the two bluer IRAC bands 
(at 3.6 /im and 4.5 /im) were merged by removing sources whose separation was smaller than 1" (roughly, 1.5 pixels 
in the mosaicked images). After the selection, we measured aperture photometry in th e 4 IRAC im a ges (fi xing the 
positions and forcing the detection in all bands) following the same technique used by IHuang et al.l (2004) . Fluxes 



were measured in small apertures of 4" diameter with SEXTRACTOR (obtaining almost identical results with other 
software, such as DAOPHOT, which was used by IHuang et al.ll2004D . The final integrated magnitude was obtained 
after applying an aperture correction based on empirical IRAC PSFs. The aperture corrections for this 4" diameter 
aperture are [0.32±0.03,0.36±0.02, 0.53±0.02,0.65±0.02] mag for channels [3.6,4.5,5.8,8.0] /im, respectively, where the 
uncertainties include the effects of typical World Coordinate System (WCS) random alignment errors (always less than 
1 pixel). For sources whose Kron aperture (in optical/NIR bands) was larger than 6" (a number chosen by studying 
our simulations described below), we measured the photometry with a large enough aperture to enclose the entire 
object and checked the results with the MAG_lSO output given by SEXTRACTOR. We would like to stress that all the 
sources in the IRAC sample have measured fluxes at both 3.6 /xm and 4.5 /im. 

The characterization of the IRAC catalogs (i.e., the analysis of the effects of confusion on the deblending of sources 
and the photometry) was carried out by simulations consisting in adding artificial sources to the mosaicked images. A 
given number of sources (7 sources/arcmin 2 , which is the number of sources corresponding to a Poisson uncertainty 
in the observed number densities) of a given magnitude were added to the images, and then the full detection and 
photometry procedure was applied. Bulge-dominated galaxies of different sizes (from 1" to 10") were also added in the 
simulation to check the photometry of nearby (z<0.5) extended objects. By measuring the angular sizes of galaxies in 
the optical/NIR images, we determined that "extended sources" (defined as sources whose semi-major axis is larger 
than 6", see below) are just a minor fraction of the total number of IRAC sources at z<0.5 (less than 3%), and 
completely negligible (less than 0.5%) at z>0.5. By checking the fraction of input sources recovered by this procedure 
(in the same position within 1" or 1.5 pixels), we estimated the completeness levels at which our catalogs are reliable 
(in terms of deblending of sources and photometry) and the accuracy of our photometry. As mentioned in Section [2j 
our IRAC catalogs are 75% complete down to [3.6]~23.3 mag. 

Our simulations also show that for sources whose semi-major axis is larger than 6", aperture photometry in a 4" 
diameter aperture corrected to an integrated flux based on empirical PSFs, underestimated (on average) the total 
flux of these sources (estimated from the MAG_iSO output given by SEXTRACTOR) in more than 10%, i.e., 1-2 times 
the typical measurement uncertainty (see below). Thus, we considered sources larger than 6" as extended sources 
in IRAC (as also recommended in the Spitzer /IRAC cookbook), and for them we estimated integrated fluxes using 
large apertures enclosing the entire objects and the extended source aperture corrections given in the Spitzer/IRAC 
cookbook. 

The errors of the IRAC photometry were estimated from the sky uncertainty (estimated with SEXTRACTOR with a 
box filtering method), detector readout noise, Poisson noise in the measured fluxes (using the detector gain and total 
exposure time per pixel), and the uncertainty in the ape rture corrections ( which include the effect of WCS errors). A 
2% absolute calibration uncertainty was also considered (jReach et al.ll2005f ). The final uncertainties were checked with 
our simulations. For each input magnitude interval, we analyzed the output magnitudes obtained with our photometric 
procedure. For [3.6]=20 mag, the typical uncertainty is 0.05 mag, and for [3.6]=24 mag, the typical uncertainty is 
0.3 mag. For [4.5]=20 mag, the typical uncertainty is 0.05 mag, and for [4.5]=24 mag, the typical uncertainty is 
0.4 mag. For [5.8]=19 mag, the typical uncertainty is 0.07 mag, and for [5.8]=23 mag, the typical uncertainty is 
0.4 mag. For [8.0]=19 mag, the typical uncertainty is 0.08 mag, and for [8.0]=22 mag, the typical uncertainty is 
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TABLE Al 

Characteristics of the data compiled for the Lockman 
Hole. 



IJcLllQ 


A cff 


m lim 


Source 








w 


IRAC-3.6 


3.561 


23.0 


Spitzer GTO 


IRAC-4.5 


4.510 


23.0 


Spitzer GTO 


IRAC-5.8 


5.689 


22.3 


Spitzer GTO 


IRAC-8.0 


7.958 


22.0 


Spitzer GTO 


MIPS-24 


23.844 


20.0 


Spitzer GTO 


B 


0.442 


26.0 


Subaru Deep imaging" 


R 


0.652 


25.4 


Subaru Deep imaging" 


I 


0.795 


25.0 


Subaru Deep imaging" 


z 


0.907 


24.5 


Subaru Deep imaging" 


U 


0.361 


23.1 


ING Archive 6 


9 


0.486 


24.0 


ING Archive 6 


i 


0.767 


22.3 


ING Archive 6 


J 


1.251 


22.5 


UKIDSS C 


H 


1.649 


20.4 


TIFKAM d 


K 


2.208 


22.9 


UKIDSS C 



Note. — (1) Name of the observing band. (2) Effective wavelength 
(in /im) of the filter+dctcctor. (3) Limiting AB magnitudes defined 
as the third quartile of the magnitude distribution of our sample. (4) 
Source from where the data were obtained: a publicly available ultra- 
deep optical data from the SMOKA Subaru Archive, taken with the 
Suprime-Cam instrument on the Subaru Telescope; b data obtained from 
the Archive of the Isaac Newton Group of Telescopes, and taken with the 
Wide Field Camera on the 2.5m Isaac Newton Telescope; c data provided 
by the UKIRT Infrared Deep Sky Survey (UKIDSS), data release 2 (DR2, 
lLawrence et al .1200711 : d data taken with the TIFKAM instrument on the 
2.1 m Telescope at Kitt Peak National Observatory. 

0.4 mag. 

All the MIPS 24 /.im data for each field (including GT O and GOODS data) were reduced and mosaicked together 
using the MIPS Data Analysis Tool ([Gordon et al.l[2005h . We detected sources and measured integrated fluxes using 
PSF fitting (with the DAOPHOT IRAF 8 package) and aperture corrections. Sources were detected in three passes to 
recover the faintest sources, many of which are hidden by brighter ones. Photometry was extracted for all the detected 
sources simultaneously. For sources of noticeable extent (more than 25"), a large enough aperture was set accordingly. 
For the rest, a circular aperture of diameter 15" (6 pixels) was utilized. For this aperture, a 17% correction in flux 
must be used to correct to the total flux (based on the theoretical PSF of MIPS). The sky estimation was carried out 
in two steps, first removing the large-scale variation (due to zodiacal light) and then measuring the background around 
each source. Based on simulations similar to those carried out with the IRAC data, we estimate that our catalogs 
are 75% completeness at F(24)=80 /xJy. Uncertainties based on these simulations are less than 5% for sources with 
F(24)>400 fiJy, and 10% for sources with F(24)~80 /xJy. 

A. 2. Optical and NIR photometry 

The Spitzer data were complemented with other publicly available and proprietary p hotometric and spectroscopi c 
data in the 3 fields. For the HDF-N and the CDF-S, the dataset is described in detail in lPerez-Gonzalez et al.l (|2005f ). 
For this paper, we added in the HDF-N the JK data described in Villar et al. (2007, in preparation; with limiting 
magnitudes 9 J=22.4 and K=21A), GALEX data extracted from th e GALEX archive (w ith limiting magnitudes 
NUV=24.9 and PW=25.3), the spectroscopic redshifts published bv iReddv etaD (|2006a[ ). and the GOODS IRAC 
and MIPS data. In the CDF-S, we added an image of size 37' x 30' taken in the NB816 filter with the Suprimc- 
Cam instrument on Subaru (with a limiting magnitude of ATB816=24.8), the spectroscopic redshifts published by 
IVanzella et ail (|2006l ). and the GOODS and Spitze r Legacy Survey (PI: van Dokkum) IRAC and MIPS data. For the 
LHF field (not used in iPerez-Gonzalez et al J 12005(1 . we summarize the main characteristics of the dataset, including 
the wavelengths, limiting magnitudes, and references for each filter in Table IAU The Subaru observations in the 
CDF-S and the LHF were obtained from the SMOKA Subaru Archive, and reduced using their pipeline SDFRED 
vl.2. The ph otometric and astrometric calibr ation was carried out by comparison with the Sloan Digital Sky Survey 
(SDSS DR4. lAdelman-McCarthv et al.ll2006l l cat alogs. The ING data were provided (fully reduced and calibrated) 
from the CASU INT W ide Field Camera Survey (jBarcons et alj|200l lYuan et alJl2003T ). The #-band TIFKAM data 
(|Le Floc'h et al.ll200l were reduced following typical NIR procedures, and the photometric and astrometric calibration 
obtained through comparison with Two Micron All Sky Survey (2MASS) catalogs. The UKIDSS data were pro vided 
(fully reduced and calibrated) by the UKIRT Infrared Deep Sky Survey (UKIDSS DR2. lLawrence et alj|2007t ). All 

8 IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research 
in Astronomy (AURA), Inc., under cooperative agreement with the National Science Foundation 

9 Defined as the third quartile of the magnitude distribution of our sample. 
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the images in our complete dataset were calibrated photometrically (using direct observations of SDSS and 2MASS 
catalogs) and astrometrically (using SDSS and 2MASS catalogs). Typical absolute photometric uncertainties were 
0.03 mag, and the WCS absolute uncertainty was always less than 0.5". 

A. 3. Merged photometric catalog 

Ape rture matched photometry in all bands was carried out using the procedure described in lPerez-Gonzalez et al.l 
( 2005). The coordinates of the IRAC detected sources are cross-correlated with each one of the UV, optical and NIR 
catalogs using a search radius of 2.5" (roughly two pixels in the original IRAC images) and starting b y the deepest 
images. Once the source was identified in one of these image (for most cases, the first one), we took the|Kron| ()1980l ) 
elliptical aperture best enclosing the entire source from this reference image, and translated it to all the other bands. 
The aperture was large enough to enclose the PSF (at least twice the FWHM of the PSF) in all UV/optical/NIR 
images (the seeing was always less than 1.5"). By randomly varying the center of this aperture in each image, we 
checked that small WCS errors did not affect the integrated apertures significantly (the variations were always well 
within the photometric uncertainties). For IRAC and MIPS, where the PSFs are comparatively large, we assumed 
the integrated magnitude measured in small apertures (applying aperture corrections), as discussed previously. For 
GALEX data, given that the FWHM of the PSF is 6" — 7" (depending on the band, position in the detector, and 
brightness of the source) we took the MAG_BEST magnitude given by SEXTRACTOR. For HST images, we picked the 
integrated flux of the closest source measured with SEXTRACTOR, not carrying out any aperture matching. For this 
reason, HST fluxes were not used in the photometric redshift and stellar mass determination. 

Uncertainties of each measured flux were obtained from the sky pixel-to-pixel variations, detector readout noise, 
Poisson noise in the measured fluxes (taking into account the detector gain and total exposure time per pixel, which 
were combined to give rms images of the fields), the errors introduced by the uncertainties in the WCS, and the 
uncertainty in the absolute photometric calibration (typically 0.03 mag). Reductions involving drizzling (e.g., in ACS 
or IRAC images), non-integer pixel shifts (e.g., NIR images), and also detector artifacts or unresolved faint sources, 
produce that uncertainties derived uniquely from pixel-to-pixel variations of adjacent sky pixels underestimate the 
real noise, since these effects correlate the signal of nearby pixels (see, e.g., ICasertano et al.l 120001 : lLabbe et al.l 120031 
iGawiser et al.|[2006tlQuadri et al.ll2007f ). To account for this, we estimated the background level and noise in 3 different 
ways. First, we measured the average signal per pixel and noise in a circular corona 5" wide surrounding the Kron 
photometric aperture for each source, scaling the noise with a N 1 / 2 factor, where N is the number of pixels of the source 
photometric aperture. To get rid of the effect of correlated noise in this estimation of the uncertainty introduced by the 
pixel-to-pixel variance, we also estimated the background level and noise using 20 artificial apertures of the same size 
as the one used for the source. These artificial apertures were built by randomly selecting (in general, non-adjacent) 
"sky pixels" in a l'xl' box around the source. Those "sky pixels" excluded the pixels whose signal was ha above 
the rms value estimated with the first method. The average signal and standard deviation of the integrated fluxes 
within these artificial apertures provided another (less biased) estimation of the background level and noise. Finally, 
we also used 20 apertures of the same size, shape, and orientation as the source photometric aperture in the 1' x 1' box 
forcing that more t han 90% of their pixels were "sky pixels" (as defined before) and applying the method described in 
lLabbe et al.l |2003) . The final background level was set to the average of the three estimations, and the background 
noise was set to the largest estimation provided by any of the three previously described methods. In practice, the 
largest estimation of the noise was, in most cases, provided by the second method: on average, the noise was 10%-20% 
higher than what was obtained with the first method (which proves the importance of correlated noise), and less than 
5% higher than the third method. 

The validity of the method used to obtain merged photometry from the UV to the MIR bands was tested by 
comparing the measured colors with those obtained from images convolved to the same resolution. For this test, we 
matched the PSFs of an optical image (the J-band) to that of the IRAC 3.6 fim channel (which is worse) using the 
IRAF PSF task (which produces a convolution kernel to match the optical PSF to the IRAC PSF). We then measured 
photometry in a 4" diameter aperture in both bands and obtain i-[3.6] colors for all the detected sources. Note that 
since both images have the same PSF, any aperture size could in principle be used to obtain colors, but very small 
apertures would be more affected by WCS and PSF matching errors. In the case of resolved nearby sources, very 
small apertures could also bias the results if the colors are not uniform across the galaxy. The colors derived with this 
method were very similar to those obtained with our photometric procedure. The absolute mean difference between 
both methods was <|A{/-[3.6]}|>=0.02 mag (the average difference was <A{i-[3.6]}>=0.005 mag), and the scatter 
0.11 mag, comparable to the color uncertainties (the average is 0.15 mag). The average difference is independent 
of the integrated magnitude and size. For sources with J<22, we find <A{/-[3.6]}>=0.004±0.12 mag, for sources 
with 22</<24 we measure <A{/-[3.6]}>=0.01±0.11 mag, and for sources with i>24 we obtain <A{7-[3.6]}>=- 
0.02±0.12 mag. For sources smaller than 6" we find <A{/-[3.6]}>=-0.01±0.12 mag, and for larger galaxies we find 
<A{/-[3.6]}>=0.006±0.11 mag. 

For some of the IRAC sources (10%-15% of the entire IRAC sample in each field), there were several UV/optical/NIR 
counterparts in ground-based images for one single IRAC source within the 2.5" search radius. For these sources, we 
remeasured the IRAC fluxes by fixing the positions of the blended objects and deconvolving the images using the 
IRAC PSFs. Although the IRAC PSFs have FWHMs of approximately 2", the determination of the central position 
of each IRAC source can be determined more accurately (the actual value depending on the brightness of the source) 
and sources are resolved for separations of the order of This means that if the source positions are known, we can 
identify and deblend IRAC sources separated ~1" from each other. We adopted a similar deconvolution method to 
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that used in lGrazian et alj (|2006f ). The ground-based optical/NIR reference image was used to measure the positions 
of the different blended sources. Then, the reference image and the IRAC images were realigned locally (in a 1' x 1' 
square region around the source) to minimize the WCS related errors in the photometry, which were expected to be 
large in the very small apertures used in the deconvolution method. The IRAC photometry in this case was measured 
by convolving the IRAC PSF with the reference image PSF and scaling the flux of each object to match the IRAC 
fluxes in an aperture of 0.9" (~1.5 pixels in the IRAC images). For this aperture size, the aperture correction of 
the IRAC bands are [1.01±0.07,1.02±0.08, 1.2±0.10,1.44±0.14] mag for channels [3.6,4.5,5.8,8.0] fxm, respectively 
(including WCS errors) . For the separations between the sources which we are trying to deconvolve (separations larger 
than 1" and smaller than 2.5"), the flux contamination from the surrounding sources to a given one was, in most cases, 
lower than a 10% of the flux in the photometric aperture. The artificial source simulations validated this procedure. 

Most of the IRAC selected sources are detected in our deepest Subaru images in the HDF-N: approximately 90% are 
detected in B, R, and/or /. In these bands, 75% of our sources are brighter than B = 25.5, R — 24.9, and / = 24.5. 
In the CDF-S, 90% of the sources are detected in the A5816 filter, and 75% of them are brighter than A5816 = 24.8. 
In the same field, about 70% of sources are detected in B (75% of them arc brighter than B = 25.3), 60% in R (75% 
of them brighter than R = 24.8), and 40% in / (75% of sources brighter than I = 23.7). MIPS at 24 fim is able to 
detect about 25% of the IRAC sources (75% of them above F(24)=40 /xJy). 

More than 90% of the /-band selected sources (see Section [5]) were also detected in deep BVRz imaging. In these 
bands, 75% of our /-band sources are brighter than B=26.0, F=25.8, #=25.5, and z=25.1. About 50%-55% of the 
/-band sample is detected by IRAC (at 3.6 and 4.5 fim; at 5.8 and 8.0 fim, the fraction drops to 40%-45%). MIPS at 
24 fxm is able to detect about 7% of the /-band sources above F(24)=80 /iJy. 

A negligible fraction of the entire IRAC sample (less than 3%, not large enough to change our results significantly) 
was detected in less than 5 filters (our limit to calculate a reliable photometric redshift), all of these galaxies presenting 
fluxes below the 75% completeness level. For the /-band selected sample, only 2% of the sources are detected in less 
than 5 filters. 

A. 4. The spectroscopic sample 

Both the HDF-N and the CDF-S include a large compilation of spectroscopic redshifts obtained by several surveys. 
Unfortunately, there is no public spectroscopi c sur vey in the LHF, so this field could not be used for building templates 
to estimate photometric redshifts (s ee Section IB .41). In t he HDF-N, we used 1,69 9 spectroscop i c redsh ifts (~20% of the 
entire sample in that field) found in lWirth et alj (|2004h . lCowie et all (|2004ft . and lReddv et alj (|2006aD . Only a fraction 
of those redshifts (1,340 sources) are flagged as high reliability (larger than 80%). In the CDF-S, we compiled 1,410 
spectroscopic redshifts (about 1 5% of the sample in that field), 891 of them flagged as reliab l e with a pr obability larger 
than 8 0%, from several sources: iLe Fevre et all ()2004l ). ls"zokolv et all (|2004l ). IVanzella et all (|2005h . and lVanzella et all 
(2006). More than half of the highly-reliable spectroscopic redshifts are below z=1.0 (55% in the CDF-S, and 80% in 
the HDF-N), and most of them are below z=1.5 (95% in the CDF-S, and 97% in the HDF-N). These spectroscopic 
redshifts were complemented with photometric redshifts estimated as explained in Section [B. 41 

A. 5. Star-galaxy separation 

In order to separate galaxies from stars in the mer ged photometric catalog s, we used eight criteria, one based 
on the STELLARITY parameter given by SEXTRACTOR (|Bertin k, Arnouti [1996) . and the other seven criteria based 
on color-color and color-magnitude diagrams using optical and NIR fluxes. All objects detected in more than one 
optical or NIR band, and presenting an average value of the STELLARITY parameter larger than 0.95 were identified 
as stars. An o bject was also con s idered a star if it satisfied any of the se col or equations ( when fluxes were available), 
extracted from lEisenhardt et all (|2004f l . iRowan-Robinson et all (|2005h . and lDaddi et alj (|2004D : a) [3.6] - [8.0] > -2 
and [3.6] - [8.0] < -1 and [8.0] < 20., or [3.6] - [4.5] > -1 and [3.6] - [4.5] < -0.5 and [4.5] < 19.5; b) [5.8] - [8.0] > -1, 
[5.8] -[4.5] < -0.2 and [8.0] < 20.; c) 7- [8.0] < -1 or/- [3.6] < 1 and [3.6] < 18. or/- [8.0] < -land [3.6] -[8.0] < -1; 
d)B-/>2x(/- [3.6]) + 0.070; e) J — K + 0.956 < 0.5; f) [3.6] 3 „ - 0.460 - [3.6] auto > -0.25 and [3.6] < 15. and 
[3.6] // — 0.460— [3. 6] auto < 0.2, or [3.6] « — [3.6] auto < —0.25, where [band\ 3 n is the magnitude in a 3" diameter aperture, 
and [band\ auto is the MAG_AUTO magnitude given by SEXTRACTOR (an estimation of the integrated magnitude); and 
g) z — A'<-0.5+0.29x(£? — z). The sta r -galaxy separation for the IRAC sample was checked a gainst the galactic 
nu mber counts published by iFazio et al.1 (|2004al see also the stellar number counts predicted bv lArendt et al.|["l998l 
and IWainscoat et"aT1ll992f ). finding very good agreement with our results (absolute differences of less than O.ldex at 
all fluxes down to the limits of our survey). Note that these authors also show that the stars dominate the number 
counts at the bright end, but they are a minor contributor at faint magnitudes (less than 4% of the sources at [3.6] >20 
are stars), the range where our extra-galactic analysis is concentrated. We have also checked that our star detection is 
able to recover more than 95% of the stars in our IRAC sample that have been spectroscopically confirmed: we identify 
222 stars out of 232 spectroscopically confirmed stars in the HDF-N, and 78 out of 82 sources in the CDF-S. All the 
objects considered stars by the spectroscopy and missed by our algorithm are extended (had effective radii larger than 
3 pixels and FWHM larger than 4 pixels) in the ACS images. In the HDF-N, our star detection algorithm identifies 6 
spectroscopically confirmed galaxies as stars, all of them being point-like in the ACS images. In the CDF-S, 14 sources 
with a spectroscopic redshift are identified as stars by our algorithm, all of them except two being point-like in the 
ACS images. 
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A. 6. AGN identification 

We used X-ray data (covering our entire surveyed regions in the 3 fields) to select candidates to harbor an AGN 
within our samp l e. In the HDF-N, we used the catalog for the 2 Ms Chandra Deep Field-North Survey published by 
lAlexander et all (|2003f) . finding an X-ray counterpart 10 for 5% of our IRAC samp le in th at field (3% of the J-band 
sample). In the CDF-S, we used the catalogs published by iGiacconi et al.l ()2002l see also lTozzi et~aT1l2006D for the 
1 Ms Chandra Deep Field-South Survey, identifying 3% of our IRAC sources as X-ray e mitters (2% of the J-b and 
sam ple). In the LHF, we ident ified AGN candidates using the XMM catalogs published by lHasinger et all (|2001l . see 
also lLehmann et al.ll2000L and lMainieri et aTll2002ah . finding an X-ray counterpart for 0.4% of the IRAC sources in 
that field (0.2% of the J-band sample). In the total IRAC sample, 3% of galaxies at any redshift were identified as 
X-ray emitters (2% of the entire J-band sample), with slightly larger values (4%-6% of all IRAC sources) found for 
sourc es at z>1.5. Observations in X-rays are known to miss very obscured AGNs (e.g., iRigbv et al.ll2004UDonlev et alJ 
2005). Other selection proced ures have been used to identify obscured AGN s, such as the presence of a power-law 
spectrum in the IRAC bands (|Alonso-Herrero et al.| [2004; Do nley et al.l l2007h . These power-law galaxies (some also 
detected in X-ray or radio data) are also a very small fraction of our IRAC selected sample, less than 1% of the total 
number of sources. We refer the reader to Section IB. 61 for a discussion about the characterization of the X-ray sources 
and the effect of AGN contamination in our results. 

B. ESTIMATION OF PHOTOMETRIC REDSHTFTS, STELLAR MASSES, AND STAR FORMATION RATES 

B.l. Stellar population synthesis models 

For the stellar population synthesis models of the SEDs of the spectroscopic sample, we carried out two sets of fits: 
1) one set assuming that the star formation history of each galaxy can be described by a declining exponential law with 
time scale r, age t (i.e., SFR(t) oc e~*/ T ), metallicity Z, and a ttenuated by an amount described by the quantity A(V) 
(1-POP models, hereafter, see also lGil de Paz fe Ma dore 2003); and 2) another set (2-POP models, hereafter) assuming 
one recent instantaneous burst of star formation of age t you , metallicity Z you and extinction A(V) you , overimposed on 
an evolved stellar population characterized by r id, toidi -^oid, and A(V) \d- The attenuation at any wavelength was 
calculated from the free parameter A(V) using the l( •harlot k- Fal (l2000l . CF00 hereafter) recipe. In this work, the 
attenuation of the gas and stellar emissions is divided into three components, based on a simple scenario: the light 
arising from the newly-formed stars, embedded in a birth cloud, is attenuated by the material in the HII region, by a 
surrounding shell of molecular and/or non-ionized atomic gas and dust, and finally by the inter-stellar medium. The 
extinction law is approximated by a power-law function of the form A\ oc A™ (the authors suggest n = —0.7). There 
is also a dependence of the birth cloud extinction on the age of the stars: for stars younger than 10 Myr (the typical 
lifetime of molecular clouds) the extinction is /x times larger than for older stars, where /z^0.3 (with significant scatter). 
We also ran a set of models assuming that the attenuation law was similar to the one found for local starbursts by 
ICalzetti et all (t2000l CALZ00 r ecipe, hereafter). The stellar em ission in our models was taken from the PEGASE 
code ()Fioc fc Rocca-Volmerangelll997l ). assuming a lSalpeterl (|1955f ) initial mass function (IMF) with 0.1<A^<100 Mq 
and a single power-law slope through the entire mass range. We also added the emission from the Hydrogen gas 
heated by the stars (emission lines and nebular continuum) using the emissio n and recombination coefficients given by 
iFerlandl ([1980) for an electron temperature T e — 10 4 K, the relations given bv iBrocklehurstJ (|1971h . and th e theoretical 
line-r atios expected for a low density gas (n c = 10 2 cm~ 3 ) with T c = 10 4 K in the recombination Case B (|Osterhrockl 
[1989T ). 

The 1-POP models required 4 parameters to fit. Our fitting routine probed the solution space in the following 
ranges for the parameters [r, t, Z, A(V)]: i) we assumed t values from an almost instantaneous burst (r = 1 Myr) to 
an almost constant SFR (r = 100 Gyr) using a logarithmic interval of O.ldex (in yr) for a total of 51 steps; ii) ages 
were probed from t — 1 Myr to t — 13.5 Gyr in logarithmic intervals for a total of 60 steps, constraining the solution 
for each object so the computed age was not larger than the age of the Universe at the redshift of the galaxy; in) we 
used the 7 discrete values of the metallicity available in the PEGASE code [0.005,0.0.02,0.2,0.4,1.0,2.5,5.0] x Z Q ; 
and iv) extinction values ranged from A(V) = mag to A(V) = 5 mag in intervals of 0.10 mag (51 steps). 

For the 2-POP models, each one of the 2 stellar populations requires in principle 4 parameters to fit, but we forced the 
recent burst to be instantaneous, so the young stellar population only requires 3 free parameters to fit. Added to those 
7 parameters, one more parameter is necessary, the burst strength 6, to describe the fraction of the total stellar mass of 
the galaxy that the recent burst has created. Our fitting routine probed the solution space in the following ranges for the 
parameters [r i d , t old , Z old , A(V) id, T y0 u, t yon , Z yon , A(V) you , b]: i) r old = 1 Myr to r old = 100 Gyr using a logarithmic 
interval of O.ldex; ii) t id — 1 Gyr to t a \d — 13.5 Gyr in logarithmic intervals (constrained by the age of the Universe at 
the redshift of each galaxy); Hi) Z old = [0.005,0.0.02,0.2,0.4, 1.0, 2.5, 5.0] xZ Q ; iv) A(V) oM = mag to A(V) i d = 5 mag 
in intervals of 0.1 mag; v) we assumed an instantaneous burst for the recent star formation event (i.e., r you = 1 Myr, 
so actually this is not a free parameter); vi) t you = 1 Myr to t you = 1 Gyr in logarithmic intervals (constrained 
by the age of the Universe at the redshift of each galaxy); vii) Z you — [0.005,0.0.02,0.2,0.4, 1.0,2.5,5.0] x Z Q ; viii) 
^4(U)oid = mag to A(V) id = 8 mag in intervals of 0.1 mag; and ix) the burst strength could take values from 0.5% 
to 15% in steps of 0.5%. 



10 Our galaxies were cross-correlated with the X-ray catalogs using a 2" search radius, as done by Rigby et al. (2004) to match sources 
at large off-axis angles in the Chandra images. 
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B.2. Stellar population synthesis fitting procedure 

The stellar population synthesis models were compared with the observed photometric data of the galaxies 
in the spectroscopic s ample using a maximum likelihood estimator similar to the one defined in Equation 6 by 
iPerez-Gonzalez et al.l (|2003bf) . which takes into account the uncertainties in each data point. All data points for 
rest-frame wavelengths bluer than 4 /jm (where stars should dominate the integrated emission of the galaxy in most 
cases) were included in the fit. 

Given the large number of possible solutions (1 x 10 6 in the 1-POP case and 3 x 10 11 for the 2-POP models), the 
amount of photometric data to fit (up to 48 bands in the case of the sources in the CDF-S, 16 in the HDF-N, and 
14 in the LHF), and the number of galaxies in our samples (more than 50,000 adding IRAC and /-band selected 
galaxies), the time requirements to probe the complete solution space for each galaxy (each one at a certain redshift) 
were prohibitively high. Therefore, we had to use a minimization procedure to search for the best solution without 
evaluating the minimization function at all poin ts in the grid of solu tions. The minimization procedure was a two step 
algorithm. First, we used a genetic algorithm ( Cha rbonneaul [19951 ) . This procedure started with 200 "individuals" 
(i.e., 200 points in the solution space), whose "genome" was formed by the 4 or 8 free parameters in our minimization 
problem. The 200 individuals were "coupled" randomly (obtaining 100 couples). Each one of these couples (formed 
by "parents" ) produced 2 "descendants" . Each descendant was built by combining randomly the parameters of the 
parents. The "genome" of the descendant had to be a better solution for the minimization problem than the "genome" 
of its parents. If, after building 10 descendants for a given couple, none or only one of them were better solutions to 
the minimization problem, the two best individuals (the best solutions to the minimization problem) were kept for the 
next generation, and the rest discarded. After every 10 combinations of parameters, we allowed a random mutation 
in one of them. After all the couples had produced 2 descendants, we eliminated the parents or descendants that 
produced the worst results for the minimization problem until 200 individuals survived, and then started again the 
procedure for another generation with the best 200 individuals. The total number of generations was set to 100. For 
the final generation, we took the 4 best individuals (the best 4 solutions of the minimization problem) and produced 
small grids of solutions around them (with a width equal to one tenth of the full size of the solution space for each free 
parameter). We evaluated all the solutions in these grids, and found the best solution and confidence intervals. Our 
minimization procedure was tested for a subsample of 1,000 galaxies in the 1-POP case by comparing the best solution 
found by the algorithm with that obtained by evaluating all the grid points in the entire solution space. For this test 
sample, the minimization algorithm recovered the best solution for ~50% of the galaxies. For the rest of sources, the 
difference between the best value and the value recovered by the minimization algorithm was always smaller or equal 
to one tenth of the size of the grid for each free parameter. We will come back to the discussion of the goodness of 
the minimization algorithm in Section IB. 51 when we discuss the quality of the derived photometric redshifts, stellar 
masses, and SFRs. 

B.3. Dust emission models 

Once the stellar spectrum was modeled, we subtracted the predicted fluxes from the photometric data points at 
rest-frame wavelengths redder than 4 jxm (if present) to obtain the emission arising from the dust. This "IR excess" 
was then fitted with one of the dust emission models of I Chary fc Elbazl (|2001). We selected the model best reproducing 
the colors of the dust emission, if several photometric points were available (for relatively low redshift; see the second 
SED fitting example in Figure [BT|) . or the model giving the closest value to the observed monochromatic luminosity if 
only one IR photometric point was available (see the third example in Figure IBlj) . To check the uncertainties in the 
derived IR-based SFRs, we also used the models of lDale fc Helo u (2002) and Rieke et al. (2007, in prepa ration) in the 
fitting of the "IR excess". The latter are empirical spect ral templates c onstr ucted largely as described in lDonlev et al.l 
(|2007f ). Appendix A. For ULIRGs, we used spectra from lArmus et a l. (2007) for the star- forming g alaxie s IRAS 1211, 
1434, 1525, 2249, in addition to the data on Arp 220 and IRAS 17208 described bv lDonlev et~aLl (|2007h . For LIRGs, 
we used IRS spectra from a MIPS GTO program led by A. Alonso-Herrero. Where the LIRGs were expected to be 
extended at the scale of the IRS slit width, they were mapped. The mapping data were reduced using "cubism" 
written by J. D. Smith as part of the SINGS legacy program 11 . The mapped spectra were collapsed into a single one 
to represent th e integrated galaxy p roperties. For all the spectra, the templates were extended to shorter wavelengths 
as described in lDonlev et alj (|2007t ); by constraining spectral segments with large beam photometry from 2MASS and 



IRAC, we are able to assemble reliable templates. Toward long wavelengths, we collected photometry from IRAS, ISO, 
Spitzer, and sub-mm facilities from NED. These data were fitted with a single blackbody with wavelength dependent 
emissivity as X 13 . 

B.4. Photometric redshifts, stellar masses, and SFRs 

Our final reference template set is composed of 2,074 galaxies (1,310 galaxies from the HDF-N and 764 from the 
CDF-S), for which we obtained 1,666 different 1-POP+dust models (each one of them with a unique combination 
of the free parameters), and 1,958 2-POP+dust models. As mentioned earlier, these galaxies were selected from the 
spectroscopic sample, all of them having a spectroscopic redshift measured with a reliability probability larger than 
80%. In addition, all the reference sources should have more than 10 different photometric data points in their S EDs 
covering the UV, optical, and NIR/MIR spectral ranges. Three examples of these templates are shown in Figure IbT] 
and discussed in Section IB. 51 The entire template set is available upon request to the authors. 



A description of the " CUBISM" software can be found at http://ssc.spitzer.caltech.edu/archanaly/contributed/cubism/ 
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The photometric redshi ft estimation for each g alaxy in our survey was carried out with our own code in a similar 
way to that described in iPerez-Gonzalez et al.l (|2005| ). Briefly, the observed data (fluxes and uncertainties) were 
compared with the red shifted models (with steps of Az=0.01) using a % 2 minimization algorithm (as the one used by 
Bolzoncll a et al.l l 2000f) . The method compared the photometry with the convolutions of the different filters with the 
redshifted templates, and determined the best template (the one giving the lowest x 2 value) for each redshift. The 
technique also included a preliminary independent detection of the 1.6 /jm bump feature (if present), which helped 
to constrain the final solution and get rid of outliers. The template giving the best solution at each redshift also had 
to provide an age of the stellar population younger than the age of the Universe at that redshift. The photometric 
redshift probability distribution was built with the best values of the x 2 estimator (corresponding to the model best 
reproducing the observed SED) at each redshift, and the most probable reds hift and uncertainty w ere estimated from 
that probability distribution (as a mean weighted with the probabilities, see iBolzonella et al.ll2000t ). 

From the best model and most probable photometric redshift, we could also obtain simultaneously an estimation of 
the stellar mass, as the model established the monochromatic luminosity per unit of stellar mass at all wavelengths. By 
scaling this model to the observed monochromatic luminosities (multiplying by a factor), we obtained the stellar mass 
of each galaxy. The final stellar mass and associated uncertainty for each galaxy were obtained as the average and 
standard deviation of the stellar masses obtained for each observed photometric band. The uncertainty includes both 
the effect of the photometric errors and the uncertainty in the determination of the redshift. These errors are estimated 
for each galaxy by considering the photometric redshift uncertainty and outliers derived from Figure [B2l for a galaxy in 
the same redshift and magnitude intervals, and studying how variations in the redshift affect the mass-to-light ratios 
in each band and the final stellar mass estimate. The average stellar mass uncertainty is 0.2dex, typical of any stellar 
population stud y (the typical accuracy of stellar masses obtained w i th stellar population synthesis models is a factor 
of 2-3; see, e.g. JPerez-Gonzalez et alj|2003d: iKauffmann et~alll2003l ; Fapovich et afll200d iFontana et al.ll2006[) . 

Star formation rates were estimated from the total IR luminosity [L(8 — 1000)] calculated by integrating the dust 
emission models for each galaxy between 8 fim and 1000 /im. Galaxies not detected by MIPS at 24 /im were assumed 
to have an upper limit flux of i^(24) = 60 /iJy. The final SF R estimation also includes the contribution from unobscured 
star formation detected directly in the UV. According to iBell et all (|2005[ ) , we can estimate the total SFR for each 
galaxy from L(8 — 1000) and L(0.28), where L(0.28) = 1^X^(0.28) is the monochromatic luminosity at 0.28 /im measu red 
d irectly from the stellar population model for each galaxy. The conversion factor is taken from iKennicuttl (|1998t ) for 
a lSalneteri <fl95fih IMF: 

SFR = 1.8 x 10~ 10 [L(8 - 1000) + 3.3L(O.28)]/L A^yr^ 1 (Bl) 

In order to characterize the uncertainties of the SFRs derived with our models, we also calculated IR-based SFRs 
by estimating monochromatic luminosities at rest-frame wavelengths 6.7 /im, 12 /im, and 15 /im. The integrated 
luminosi ty L(8 — 1000) can be o btained from these monochromatic luminosities by applying the empirical relationships 
found in lCharv fc E lbaz ( 20011). Another estimatio n of the SFR can be obtained from the rest-frame monochromatic 
luminosity at 24 /im applying the equation given in Alon so-Herrero et a l. (2006). We will discuss the uncertainties in 
the SFR estimations in Section lB.5.41 

B.5. Evaluation of the modelling procedure and derived parameters 
B.5.1. Some examples of SED fits 

Figure [Bl] shows three examples of the dust and stellar population models for IRAC sources at different redshifts. 
The three panels on the left show the fits for the 1-POP case, and the three panels on the right show the fits for the 
2-POP case for the same sources. 

The upper two panels present a source nicely fitted by a single old stellar population with intermediate extinction, no 
current star formation (a bulge dominated galaxy), and no detection at 24 /im. In this case, although the photometry 
at rest-frame wavelengths redder than 4 /im was not used in the model fitting, those points are well reproduced with 
just stellar emission. When fitting the same SED with a 2-POP model, we recover very similar parameters to the 
1-POP case, with a minor contribution (just 1% in mass) from a more recent burst. Note that both types of models 
give very similar stellar mass values. 

The second example (middle row) shows an intermediate redshift galaxy detected at 24 /mi. This galaxy can be fitted 
either by an old single stellar population with a large extinction or with a combination of an old stellar population 
with very low attenuation and a more recent burst contributing about 13% to the total stellar mass of the galaxy. This 
recent burst presents a relatively large dust attenuation that could be responsible for the emission in the MIR/FIR 
(the galaxy is detected at 24 /im). For this galaxy, dust emits a significant fraction (about 50%) of the 8 /im luminosity 
(rest-frame 4 /im) and almost 100% of the 24 /im luminosity (rest-frame 12 /im). The 1-POP models give a larger 
stellar mass value than the 2-POP models (still, the difference is a factor of 3, comparable to the typical uncertainty 
in stellar population studies) because a lot of stars are necessary to fit the high NIR photometric data points, while a 
lot of extinction is necessary to simultaneously fit the UV/optical fluxes. 

The third example (bottom row) is a high redshift galaxy with a very blue spectrum. It can be fitted either with an 
almost continuous unattenuated star formation (based on the high r value) lasting about 100 Myr, or with a similar 
primary burst (producing 93% of all the stellar mass) followed by a more recent (10 Myr) and much more attenuated 
(A(V) you —7 mag and a strong MIR emission detected at 24 /im) event of star formation. Note also that the IRAC 
photometric point at 8.0 /im is too high in comparison with the combined star+dust models. At wavelengths around 
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Fig. Bl. — Three examples of the stellar population and dust emission modelling of IRAC selected sources in the spectroscopic sample. 
The spectroscopic redshift and main stellar population parameters of the best fit are given in each panel. Filled black stars and vertical 
error bars show the photometric points used in the stellar population fits (wavelengths bluer than 4 fim). Horizontal error bars for each 
photometric point show the width of the filter. Open black stars are the photometric data points used in the modelling of the dust emission. 
The left panel of each row shows the 1-POP stellar emission fits with a cyan line, and the final fit (including nebular continuum and emission 
lines) with a red continuous line. On the right panel of each row, the same SEDs have been fitted with the 2-POP models, where one stellar 
population is plotted with a cyan line, the other population with a magenta line, an d the addition of both with a red line (including nebular 
continuum and emission lines). For all panels, the dust emission model taken from Chary & Elbaz (2001) which best reproduces the MIR 
emission (if present) has been plotted with a dashed red line. Green vertical lines show the positions of the most interesting emission-lines 
in the optical and NIR spectral ranges. 
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Fig. B2. — Comparison of the spectroscopic and photometric redshifts for IRAC selected sources in the HDF-N (top panel) and the 
CDF-S (bottom panel). Gray symbols are sources with spectroscopic redshifts which have a reliability probability lower than 80%. Open 
stars are sources detected in less than five bands. The dashed lines show the equality line, and the <r z /(l+z)<0.1 area. 



A~4-10 /im or even at A~2-10 /mi for very luminous IR sources (with very hot dust), the integrated emission comes 
from both the dust and the stars in comparable fractions. In this overlap region between the dust and stellar models, 
the spectrum may show prominent emission-lines, PAH features, or emission from hot dust (e.g., coming from a dust 
torus surrounding a nuclear massive black hole) which are not found in the stellar and du st emission models. F or 
example, there is a PAH feature at rest-frame 3.3 /mi (very weak in all dust models in the ICharv fc Elbaz 2 0011 or 
Dale & Helou 2002 libraries) which may have a non-negligible contribution to the global emission in this spectral region 
for very luminous IR sources. 
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B.5.2. Statistical evaluation of the photometric redshifts 

The main three parameters that we want to extract from the SED fits are the photometric redshift, the stellar mass, 
and the SFR of each galaxy. The quality of our photometric redshifts is checked in Figure IB2I for the fields with 
extensive spectroscopic data: the HDF-N and the CDF-S. Unfortunately, given that there is not a systematic public 
spectroscopic survey in the LHF, we cannot check our photometric redshifts directly in this field. In spite of this, the 
photometry in the LHF is as good or even better (given that the optical images are ultra-deep observations taken with 
Subaru) than in the other two fields, and the general redshift distribution for the LHF sources is similar to that in the 
HDF-N and the CDF-S. Consequently, we conclude that the quality of the photometric redshifts in the LHF must be 
comparable to the other two fields (see also the discussion about the redshift distribution of our sample in Section 0]). 

The top panel of Figure IB2I shows the comparison of our photometric redshifts and spectroscopic redshifts for 
the IRAC selected sample in the HDF-N (for the 1,702 sources with available spectroscopy). The average (median) 
redshift difference (<5z=z spcc -z p hoto) is 0.014 (0.010), comparable to the redshift step used in our photometric redshift 
technique. This demonstrates that there are no systematic errors in our redshifts. Almost all sources, 95%, have values 
of <7 z /(l+z)<0.2 (where er z is the absolute value of Sz), 88% of the objects present values of <7 z /(l+z)<0.1, and 70% 
have <7 z /(l+z)<0.05. The average (median) value of <r z /(l+z) is 0.055 (0.032). Very similar statistics are obtained for 
the /-band selected sample: 94% of these sources present <r z /(l+z)<0.2, 86% cr z /(l+z)<0.1, and the average (median) 
cr z /(l+z) is 0.060 (0.036). 

The quality of the photometric redshifts in the CDF- S for IRAC selected sources (for the 1,410 sources with available 
spectroscopy) is shown in the bottom panel of Figure lB2l The average (median) value 5z is 0.020 (0.015). In this 
field, 93% of the objects have values of <r z /(l+z)<0.2, 85% of the objects have values of cr z /(l+z)<0.1, and 67% have 
£r z /(l+z)<0.05. The average (median) <r z /(l+z) is 0.060 (0.040). For the /-band selected sources, the numbers are 
similar: 92% of these sources present er z /(l+z)<0.2, 80% ct z /(1+z)<0.1, and the average (median) ct z /(1+z) is 0.080 
(0.047). 

Some of the sources used in the photometric redshift evaluation depicted in Figure IB2I were used in the building 
of the templates (77% of all sources plotted in the HDF-N, and 54% in the CDF-S). The validity of our procedure 
(including the merging of photometric data) and the templates built with data from the HDF-N and the CDF-S was 
also tested in completely different and independent fields. For example, in the Extended Groth Strip (Perez-Gonzalez 
et al. 2007, in preparation), we compared our photometric redshifts (based on photometry measured in the same way 
as in this paper) with spectroscopic values for 6,828 sources, obtaining that for 87% of the galaxies, our photometric 
redshifts wer e better than cr z /(l+z)<0.1, and for 95% were better than er z /(l+z)<0.2. 

Figure lB2l demonstrates the high quality of our photometric redshifts at z<1.5. Beyond this redshift, spectroscopic 
surveys have severe limitations due to the intrinsic faintness of the sources (most of them are below the typical R^25 
spectroscopic limit) and the absence of bright spectroscopic features in the observed optical range for sources at 
1.5<z<2.5 (the redshift desert). Therefore, photometric redshifts cannot be extensively tested at high-z, given that 
very few spectroscopic redshifts are available. To overcome this problem as much as possible, we included up to 59 
galaxies at z>1.5 in our template set, most o f them extracted from spectroscopic surv eys carried out with spectrographs 
with enhanced sensitivity in the blue (e.g.. ISteidel et al1 l2b04: Re ddv et al.ll2006al ). Using the very few sources with 
reliable spectroscopic redshifts at z>1.5, our photometric redshifts seem to degrade to some extent. In the HDF-N, 
69% of the sources at z>1.5 have photometric redshifts er z /(l+z)<0.2 and 50% with u z /(l+z)<0.1. In the CDF-S, 
86% of the sources at z>1.5 have photometric redshifts <7 z /(l+z)<0.2 and 59% with <t z /(1+z)<0.1. 

These statistics are highly biased against red objects, and for blue sources, a very small number of sources is used in 
the comparison (less than 30 in each field). To further test our results at high-z, we analyzed the photometric redshift 
distribution of samples of galaxies selected with the different color techniques described in Section [6] We considered all 
the IRAC sources identified as LBGs with i?<25.5, and DRGs or/and BzK galaxies with if <22.9 (if[Vega]<21). In 
Section [SI we demonstrate that our IRAC survey detects virtually all these sources (given that the number densities in 
our sample are very similar to those measured by other surveys focused on the detection of these high-z populations) . 
Here, we test the photometric redshifts derived for these sources, a topic that will be discussed in more detail in a 
future paper (Barro et al., 2007, in preparation). 

The average redshift of the LBG-BM sources in our IRAC sample is <z>=1.4±0.3. Both the average and standard 
deviation values agree, within the typical p hot ometric redsh i ft unce rtainties, with the average spectroscopic value of 
<z>=1.7±0.3 given bv ISteidel et~a!] (|2004h . In lReddv et"all (|2006al ). they also find an average <z>=1.7±0.3 for the 
LBG-BM sources in the HDF-N, where our average photometric redshift is <z>=1.6±0.3. If we calculate the average 
spectroscopic redshift for the LBG-BM sources in our IRAC sample with available spectroscopy (9% of the total), we 
obtain <z>=1.3±0.3. Our average is also consistent with the average photometric redshift published bv lQuadri et alJ 
(|2007f ) for LBG-BM galaxies in MUSYC, <z>=1.4, a nd the first peak of the photometric redshift distribution of LBGs 
in the GOODS-MUSIC sample (|Grazian et al.|[2007h . also at z~1.4. 

The average redshift of the LB G-BX sources in our IRAC sample is <z> =2.0±0.4, consistent with the spectroscopic 
values of <z>=2.2±0.3 found bv lSteidel et~aT1 (|2004f ) and <z>=2.2±0.4 by lReddv etldl (|2006al ) in the HDF-N (where 
we obtain <z>=2.1±0.3). The LBG-BX sources in our sample with available spectroscopy (5% of this sub-sample) 
have an average spectroscopic reds hift of <z>=1 . 7±0.4 . Our average is again in perfect agreement with the average 
photometric redshift published bv iQuadri et al.l (l2007f) for this p opulation, <z>=2.1, and the second peak of the 
photometric redshift distribution of LBGs in lGrazian et al.l (|2007f ). placed at z^2.2. 

The sources in our sample identified as "classical" LBGs lie at an average redshift of <z>=3.1±0.5, which compares 
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well with the spectroscopic values from ISteidel et all (|2003l) . ISteidel et all (|2004f ) and iReddv et all (|2006af) . all of 
them being <z>=3.0±0.3. Only 2% of our sources identified as "classical" LBGs have spectroscopy, and the average 
spectroscopic redshift for them is <z>=2.5±1.0. 

The average photometric redshift of the population of DRGs identified in our IRAC survey is <z>=2.2±1.0, and the 
median is z=2.5, in good agreement (taking into account photometric redshift uncertainties) with the average spec- 
trosco pic redshift <z>=2.5±0.4 in the HDF-N (we obtain <z>=2.4±0. 9 just in this field) published bv IReddv et al.l 
(|2005l) . the median and rms photom e tric v alues z=2.6±0.7 published bv iFranx et al.l (|2003[ ), th e median photometri c 
redshift z=2.2 from iPapovich et alj (|2006t ). and the median photometric redshift z~2.5 from iQuadri et aT] (|2007f ). 
Spectroscopic redshifts are available for just 4% of the DRGs in our sample, with an average of <z> = 1.5±0.9, a 
lower value than the photometric estimation, probably due to the bias of spectroscopic surveys towards the optically 
brightest sources (whose probability of being at lower redshifts is relatively larger). 

BzK (combining both PE and SF sub-types) sources in our IRAC sample have an average photometric r edshif t 
<z>=2.1±0.6, which compares nicely with the average spectroscopic value <z>=2.1±0.4 fromlReddy et al.l (l2005h. 
Other photometric redshift studies obtain similar redshift distribution, e.g.. IQuadri et all (|2007[ ) and iGrazian et al] 
(2007). The average spectroscopic redshift for BzK sources in our sample (available just for a 1% of the total number 
of BzK galaxies) is <z> = 1.7±0.3. 

The previous statistics and the consistency with spectroscopic and photometric values found in the literature demon- 
strate that our photometric redshifts for the galaxies at z>1.5 are also reliable. Still, a spectroscopic survey focused 
on high redshift sources is necessary to increase the reliability (narrow the uncertainties) of our results at z>1.5. 

B.5.3. Statistical evaluation of the stellar masses 

In this section, we discuss the quality of our stellar masses, and the possible systematics introduced by our fitting 
algorithm and the a priori assumptions of the models. 

First, we checked how well the minimization algorithm of our SED fitting technique recovered the stellar mass value 
corresponding to the model best fitting the data. For that purpose, we used 1,000 randomly selected galaxies for 
which we probed all the nodes in the solution grid for the 1-POP case. On average, the difference between the stellar 
mass estimated with the minimization algorithm and the stellar mass given by the model best fitting the SED is 
0.002dex, the median is O.OOOdex, the standard deviation is 0.07dex, and there are not any absolute differences larger 
than 0.20dex. For the 2-POP case, the number of points in the solution grid is too large to attempt the individual 
evaluation of all of them. To test this case, we only considered 100 randomly selected galaxies and degraded the 
resolution of the parameter space grid by one third for all the free parameters (thus, we only considered 1 x 10 8 
models). On average, the difference between the stellar mass estimated with the minimization algorithm and the 
stellar mass given by the 2-POP model (with a coarse solution grid) best fitting the SED is 0.04dex, the median is 
0.02dex, the standard deviation is 0.15dex, and there are not any absolute differences larger than 0.30dex. These 
statistics confirm that the minimization algorithm is able to recover the best stellar mass estimate within the typical 
uncertainties in stellar populations synthesis analysis (a factor of 2-3). 

We also compared the stellar masses obtained with the 1-POP and 2-POP models. For about 70% (55%) of the 
galaxies, both estimates are equal within a factor of 0.3dex (0.2dex). However, for the rest of galaxies (virtually all of 
them with A4<10 105 A4q), the 2-POP estimates are higher (with the most extreme cases showing a difference of up to 
a factor of 10). On average, including all galaxies, stellar masses derived with 2-POP models are 0.18dex higher than 
those derived with 1-POP models. For galaxies with A / l>10 10 ' 5 A4q, the average difference is significantly smaller, 
just 0.02dex (with a scatter of 0.15dex). This can be explained by the fact that most of the photometric data points 
in the modelling fits are found in UV/optical wavelengths, where the emission of relatively young stars is significant. 
Older stars, possibly much more numerous and dominating the global stellar mass of a galaxy, may be hidden by the 
intensity of more recent starbursts. This effect should be more noticeable in less massive systems presenting bright 
recent bursts involving a relatively high fraction (larger than what is normally observed in very massive galaxies with 
old stellar populations) of the total stellar mass of the galaxy. Only in the 2-POP models are we able to take this 
effect into account, and that is why in this case we systematically obtain larger stellar masses for some galaxies with 
relatively low masses. We conclude that the choice of the 1-POP or 2-POP models does not change the stellar masses 
significantly (more than the typical uncertainties of a factor of 2-3) in a statistical sense, and the effect on the masses 
for massive galaxies (which dominate the stellar mass density at any redshift) is very small. 

Uncertainties in the stellar emi ssion models are known t o introduce systematic errors in the estimation of stellar 
masses from photometry (see, e.g. .Ivan der Wei et alJl2006allbT) . In order to check this effect, the stellar masses obtained 
with the PE GASE code dFioc &: Rocca-Vormerangelll997l ) were compared with the values estimated by using the BC03 
models from lBruzual fc Chariot) (|2003h . O n average, the BC03 models give stellar masses larger by 0.03dex (less than 
10%), with a scatter of 0.18dex. For 95% of the galaxi es, the ste l lar m ass diffe rence is lower than a factor of 3. We 
also fitted the SEDs with the M05 models developed bv lMarastonl (|2005l see also lBruz ual 2007)), which include a more 
refined treatment of the emission from thermally pulsating asymptotic giant branch stars, and are claimed to obtain 
stellar masses that can be lower by as much as 60% (based on the prediction of lower NIR mass-to-light ratios for some 
ages). On average, the M05 models give stellar masses smaller by 0.14dex (less than 30%), with a scatter of 0.22dex 
(and no clear dependence on redshift). For 96% of the galaxies, the stellar mass difference is lower than a factor of 3. 

One important a priori assumption of any stellar population modelling is the treatment of extinction by dust. We 
compared the stellar masses obtained with the two different extinction recipes we considered (CF00 and CALZ00). 
For about 80% (65%) of the galaxies, both estimates are equal within a factor of 0.3dex (0.2dex). For the rest of 
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galaxies (again, most of them with .M<10 10 ' 5 AAq), the estimates using the CFOO recipe are higher up to a factor of 5. 
On average, stellar masses derived with th e CFOO recipe are O.lOdex higher than those derived with CALZOO law. As 
discussed in lPerez-G onzalc z~et al.l (|2003cf) , in the CFOO recipe the attenuation of the emission arising from the stars 
is always (except for very young bursts) larger than the attenuation of the gas emission. The CALZOO recipe shows 
a opposite behavior, given that the attenuation of the stellar emission is roughly half of the attenuation of the gas 
emission. Moreover, the attenuation wavelength dependence (from the UV to the NIR) proposed by CFOO is shallower 
than the one in CALZOO. This leads to a need of more stars to obtain the same observed luminosity for equal values of 
the extinction in the CFOO case, which explains the larger stellar masses derived for this case (on average). However, 
the final effect on the masses is of the order of O.ldex, which demonstrates that choice of an extinction recipe does not 
change the stellar masses more than the typical uncertainties. 

Another important assumption of the stellar population models is the IMF, which has a direct effect on the derived 
stellar masses. Different IMFs produce stellar spectra with very similar colors, b ut with less or more stars, which 
causes a sy stematic uncertain ty in the final stellar mass estimations. For example, a lKroupa et ab | (1993) I MF (as the 
one us ed in lBorch et al.|[2006l ) pred icts stellar masses small er than ours by a factor of ~1.7, or a lBaldrv k, Glazebrookl 
(2003) IMF (used in, for example. [Glazebro ok et aTll2004h predict s also smaller m asses by a factor of ~1.8. All our 
results and those extracted from the literature were normalized to a lSalpeterl (1955) IMF with 0.1<A / (<100 A4q. If the 
IMF is universal (the same at all redshifts), this choice should not affect our results other than an overall normalization. 
A discussion about changes in the IMF from galaxy to galaxy is far out of scope of this paper. 

Finally, we performed another test of the goodness of our stellar mass estimates by comparing the results obtained 
from direct comparison of the SEDs with the entire grid of stellar population models (once the redshift of a galaxy is 
known) with the results obtained with the photometric redshift technique using the empirically built set of models, 
from which we obtained stellar mass estimates for all galaxies. We find a very good agreement between these two 
stellar mass calculations: 90% of galaxies present an average difference between the two mass estimates of less than 
O.Oldex, and the scatter around this value is 0.15dex. 

Based on this discussion, the choices of 1-POP or 2-POP models, distinct stellar population libraries, different 
IMFs, or different extinction recipes produce changes in the derived stellar masses of the same order or smaller than 
the typical error in any stellar population synthesis analysis (a factor of 2-3), directly linked to the degeneracies of 
the solutions to the problem. Thus, in the Se ctions [H to [71 we will only present the resul ts obtained with the stellar 
masses estimated with the 1-POP models, the ICalzetti et al.l (j2000t ) extinction law, and a ISalpeteil (|1955h IMF. This 
choice will also allow us to compare directly with other previous works found in the literature, that usually assume 
these characteristics in their modelling procedures. 

B.5.4. Statistical evaluation of the SFRs 

In order to understand the systematic and random uncertainties of our estimations of the SFR for each galaxy, we 
carried out two tests. 

First, we used 3 different dust emission template sets built by ICharv fc Elbazl (|200lD . iDale k Heloul (|2002f ). and 
Rieke et al. (2007, in pre paration). The va lues of the IR SFR [est imated from L(8 — 1000) using the conversion factor 
found in iKennicutt 1998] d eriv e d wit h the ICharv fc E lbaz (2001) models were systematically smaller than the SFRs 
derived with the IDale fc Heloul (|2002f ) models (on average, O.ldex) and Rieke et al. (2007, in preparation) templates 
(on average, 0.2dex). To take into account the systematic uncertainties introduced by the use of a particular set of 
models, we finally considered an average value of the estimations from the three template sets. The typical uncertainty 
of this average value (based on the standard deviation of the 3 estimations) is about 50%. 

The second test consisted in obtaining IR-based SFRs with different methods. Classically, IR-based SFRs are 
calculated from the integrated IR luminosity L(8 — 1000). The quantity L(8 — 1000) can be estimated for each 
galaxy by fitting the IR spectrum with models of the dust emission. For our galaxies, this translates to a significant 
extrapolation in the templates, since the reddest point in our SEDs corresponds to the observed MIPS 24 fim emission, 
and we are assuming that a color or a single photometric point in the MIR is closely related to the emission in the 
FIR, which dominates the integrated IR luminosity. However, one can also avoid this large extrapolation by estimating 
monochromatic luminosities at specific wavelengths which are not far from the reddest photometric point in our SEDs. 
In this sense, we estimated monochromatic luminosities at 6.7 /im, 12 £tm, and 15 iiyh, and the n calculated the 
integrated luminosities L(8 — 1000) using the empirical relationships built by ICharv &: Elbazl ((20011 ) . Since they are 
based in the same templates, these estimations of L(8 — 1000) based on different monochromatic emissions are not 
independent. However, another independent SFR estimation was obtained by extrapolating in the models to measure 
the rest- frame MIPS 24 ^m monoch romatic luminosity. This luminosity was converted to a SFR using the calibration 
given in lAlonso-Herrero et a l. (2006). The typical scatter of these different IR-based SFR estimations obtained from 
monochromatic emissions is 30%. 

From these tests, we conclude that our IR-ba sed SFR estimations a r e good within a facto r of 2, which is cons istent 
with other evaluations of IR-based SFRs (e.g., iPapovich fc Bel]||2002t ILe Floc'h et"aIH2005t ICaputi et al.ll2006l) . 

B.6. Evaluation of parameters derived for galaxies with AGNs 

According to Section IA.61 a small fraction (less than 5%) of our galaxy sample probably harbor an AGN which 
emits strongly at X-ray and/or IR wavelengths. Given that our goal is to estimate the total stellar mass content of 
the Universe at any redshift up to z~4, we must try to keep this type of sources in our sample. However, the emission 
of the dust heated by the nuclear massive black hole can extend to the NIR (if very hot dust, with a temperature of 
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T>100 K, is present) and even to optical bands (e.g., in the case of Type 1 QSOs), affecting the estimations of the 
photometric redshifts, stellar masses, and SFRs. 

Given that our photometric redshifts are mainly based on stellar population synthesis models, we can expect that 
galaxies whose UV-to-NIR SEDs are not stellar are not well represented by our template set, and there is a large 
probability that the photometric redshift estimation fails. However, only the most extreme and powerful AGNs in our 
sample would affect the UV-to-NIR global SED of the host galaxy. This is demonstrated bv lDonlev et al.l ((2007), who 
build median rest-frame SEDs of X-ray-detected IRAC sources in the HDF-N, finding that only galaxies with X-ray 
observed luminosities (integrated from 0.5 to 8 keV) L(X)>10 44 erg s _1 present non-stellar SEDs, and galaxies with 
L(X)=10 43-44 erg s _1 start to show significant emission from hot dust at A rest >2 /iin. Consequently, our stellar mass 
estimates should only be affected by the presence of an AGN for bright X-ray sources. To further test our stellar 
masses for AGN, we run a set of stellar population models on 1,000 randomly selected galaxies fitting the SED only 
up to A rcs t=2.3 (im (the if-band) instead of up to A rcs t=4 fim, to exclude the hot dust emission that can arise at 
A rcst ~2-4 /xm (note that only one photometry data point at most is removed from our SEDs in this case). The average 
difference between the masses estimated in this way and our original values is just 0.002dex, and the scatter is O.lOdex. 
Cutting the SEDs at even bluer wavelengths, A rcs t = 1.5 (im has a similar negligible effect: the average difference is 
0.004dex, and the scatter is 0.13dex. This demonstrates that the (possible) AGN emission in the NIR does not bias our 
stellar mass estimates for X-ray sources of moderate brightness. For the brightest sources, however, the UV-to-MIR 
SED is significantly affected by the AGN emission, so we decided to remove from our sample all the X-ray detected 
sources with L(AT)>10 44 erg s _1 . These sources are just 0.4% of the entire IRAC sample (0.5% of the sample in the 
HDF-N and the CDF-S, and 0.1% in the LHF), slightly more common at z>2 (l%-2% of all IRAC selected galaxies 
at high redshift), so they should have a very small additive contribution (not accounted in our results) on the stellar 
mass functions and densities. 

Concerning the estimation of photometric redshifts, their reliability for X-ray sources is slightly lower than for the 
global sample: we obtain redshifts with er z /(l+z)<0.1 for 81% of the X-ray detections, and with cr z /(l+z)<0.2 for 
92%. As mentioned earlier, these sources are very few in comparison with the entire IRAC sample, and they are 
accounted in the stellar mass function estimate by using the photometric redshift uncertainties (which are estimated 
without removing this type of sources). 

In the case of the SFR estimations, dust obscured AGNs are expected to radiate the absorbed emission in the 
MIR/FIR. Although it is common that star formation coexists with AGN activity, it is not possible to decompose the 
IR emission into the components coming from the two different phenomena, mostly when very few photometric points 
are available in this wavelength range. Therefore, all X-ray sources were removed in the analysis of the (specific) SFRs 
carried out in Section [7.21 Given that we were only interested in the distribution of specific SFRs of our sample, our 
results are not significantly affected by the exclusion of this type of sources. 



